

//注意到，访问矩阵内部的元素，如m.data[i][j],其实相当于一个三级指针。
//相当与用m的指针，加上偏移量找到data的位置，再加上偏移找到i行的位置，最后偏移到j列。中间有3次偏移过程。
//如果定义real_t **m_data=m.data;用m_[i][j]访问数据，则访问数据可以提高1/3的速度，这里的优化是很可观的。

//现在代码已经写完了，如果人工在每个函数前给每个矩阵定义一个二级指针以提高访问速度，效率太低。
//应该在编译器上做优化。可以先做文本上的优化。
//比如：检查每个函数，记录每个多级指针的调用频率，如果调用频率比较高，就设置新指针。



#include<math.h>
#include<stdlib.h>

#include"sl_matrix.h"




//矢量内存分配。 v：矢量对象（这里必须用引用，因为要修改指针指向）。size：矢量大小。
void vec_allocate(vec* v, int size)
{
	v->size = size;
	v->data = (real_t*)malloc(size * sizeof(real_t));

}
//回收矢量内存
void vec_deallocate(vec* v) 
{
	free(v->data);
	v->data = NULL;
}


//矢量内存分配。 v：矢量对象（这里必须用引用，因为要修改指针指向）。size：矢量大小。
void ivec_allocate(ivec* v, int size)
{
	v->size = size;
	v->data = (int*)malloc(size * sizeof(int));

}
//回收矢量内存
void ivec_deallocate(ivec* v)
{
	free(v->data);
	v->data = NULL;
}


//内存分配。 
void mat_allocate(mat* m, int rows, int cols)
{
	m->cols = cols;
	m->rows = rows;
	m->data= (real_t**)malloc(rows * sizeof(real_t*));
	for (int i = 0; i < rows; i++)
	{
		m->data[i] = (real_t*)malloc(cols * sizeof(real_t));
	}
}
//回收内存
void mat_deallocate(mat* m)
{
	for (int i = 0; i < m->rows; i++)
	{
		free(m->data[i]);
	}
	free(m->data);
	m->data = NULL;
}
//

/// <summary>
/// 向量打印函数
/// </summary>
/// <param name="v">打印的向量</param>
void vec_print(vec v )
{
	real_t* data = v.data;
	int i;
	printfs(" \r\n vec: [ ");
	for (i = 0; i < v.size; i++)
	{
		printfs(" %9f , ", data[i]);
	}
	printfs(" ]  ");

}

/// <summary>
/// 向量打印函数
/// </summary>
/// <param name="v">打印的向量</param>
void ivec_print(ivec v)
{
	int* data = v.data;
	int i;
	printfs(" \r\n ivec: [ ");
	for (i = 0; i < v.size; i++)
	{
		printfs(" %9d , ", data[i]);
	}
	printfs(" ]  ");

}
/// <summary>
/// 控制台打印矩阵
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="precision">数据保留有效数字</param>
/// <param name="data_width">每个数字的宽度</param>
/// <returns>true</returns>
void mat_print(mat m)
{

	int i;
	printfs(" matrix: ");
	for (i = 0; i < m.rows; i++)
	{
		printfs(" \r\n [ ");
		for (int j = 0; j < m.cols; j++) {
			printfs(" %9f , ", m.data[i][j]);

		}
		printfs(" ]  ");
	}
	printfs(" \r\n  ");


}


/// <summary>
/// 向量随机化
/// </summary>
/// <param name="v">待随机化的向量</param>
/// <param name="seed">随机数种子</param>
/// <param name="min">随机数的最小值</param>
/// <param name="max">随机数的最大值</param>
void vec_randomize(vec v, int seed, real_t min , real_t max )
{
	real_t* data = v.data;

	int meisen = 2147483647;   //m是梅森素数（形式为2^(p)-1的素数）：m=2^(31)-1  =  2147483647
	int a7_5 = 16807;		  //a=7^(5)
	unsigned long int x_i = seed;  //0-m之间的整数。是余数。
	real_t random_number;

	x_i = (a7_5 * x_i) % meisen;
	random_number = (real_t)x_i / meisen;  //现在random_number是0-1的分布	
	random_number = min + (max - min) * random_number; //把random_number缩放到min到max的分布


	for (int i = 0; i < v.size; i++) { //用随机数初始化向量
		x_i = (a7_5 * x_i) % meisen;
		random_number = (real_t)x_i / meisen;  //现在random_number是0-1的分布	
		random_number = min + (max - min) * random_number; //把random_number缩放到min到max的分布
		data[i] = random_number;
	}
}
	

/// <summary>
/// 向量模长归一化
/// </summary>
/// <param name="v">归一化的向量</param>
void vec_normalize(vec v)
{
	real_t* data = v.data;
	real_t f = 0;
	for (int i = 0; i < v.size; i++) {
		f += data[i] * data[i];
	}
	f = sqrtf(f);
	for (int i = 0; i < v.size; i++) {
		data[i] = data[i] / f;
	}
}

/// <summary>
/// 向量初始化
/// </summary>
/// <param name="v">初始化的向量</param>
/// <param name="r">初始化的值</param>
void vec_init(vec v, real_t r )
{
	int i;
	real_t* data = v.data;
	for (i = 0; i < v.size; i++)
	{
		data[i] = r;
	}
}

/// <summary>
/// 实现向量的复制。 A >> B
/// </summary>
/// <param name="vA">被复制的向量</param>
/// <param name="vB">待复制的向量</param>
/// <returns>若vA和vB是等长的向量，可以复制，返回true。否则返回false</returns>
bool vec_copyto(vec vA, vec vB)
{
	real_t* dataA = vA.data;
	real_t* dataB = vB.data;
	if (vA.size == vB.size)
	{
		for (int i = 0; i < vA.size; i++) {
			dataB[i] = dataA[i];
		}
		return true;
	}
	else return false;
}

/// <summary>
/// 实现向量的复制A << B
/// </summary>
/// <param name="vA">待复制(赋值)的向量</param>
/// <param name="vB">被复制的向量</param>
/// <returns>若vA和vB是等长的向量，可以复制，返回true。否则返回false</returns>
bool vec_copyfrom(vec vA, vec vB)
{
	real_t* dataA = vA.data;
	real_t* dataB = vB.data;
	if (vA.size == vB.size)
	{
		for (int i = 0; i < vA.size; i++) {
			dataA[i] = dataB[i];
		}
		return true;
	}
	else return false;
}

/// <summary>
/// 判断两个向量是否相等
/// </summary>
/// <param name="vA">向量A</param>
/// <param name="vB">向量B</param>
/// <param name="precision">精度要求</param>
/// <returns>如果向量每个元素的差都小于精度要求，说明两个向量相等</returns>
bool vec_isEqual(vec vA, vec vB, real_t precision )
{
	real_t* dataA = vA.data;
	real_t* dataB = vB.data;
	if (vB.size == vA.size)
	{
		for (int i = 0; i < vA.size; i++) {
			if (ABS(dataA[i] - dataB[i]) > precision) return false;
		}
		return true;
	}
	else return false;
}

/// <summary>
/// 判断向量是否零向量
/// </summary>
/// <param name="v">向量</param>
/// <param name="precision">精度要求</param>
/// <returns>如果向量每个元素的值都小于精度要求，说明是零向量</returns>
bool vec_isZero(vec v, real_t precision )
{
	real_t* data = v.data;
	for (int i = 0; i < v.size; i++) {
		if (ABS(data[i]) > precision) return false;
	}
	return true;
}

/// <summary>
/// 实现向量之间的加法
/// </summary>
/// <param name="A">向量A</param>
/// <param name="B">向量B</param>
/// <param name="C">向量C=A+B</param>
/// <returns>运算成功返回true。如果向量ABC的元素个数不一致，运算失败，返回false。</returns>
bool vec_add(vec A, vec B, vec C)
{
	real_t* A_data = A.data;
	real_t* B_data = B.data;
	real_t* C_data = C.data;
	if ((A.size == B.size) && (A.size == C.size))
	{
		for (int i = 0; i < A.size; i++) {
			C_data[i] = (A_data[i] + B_data[i]);
		}
		return true;
	}
	else {
		return false;
	}
}

/// <summary>
/// 实现向量之间的减法
/// </summary>
/// <param name="A">向量A</param>
/// <param name="B">向量B</param>
/// <param name="C">向量C=A-B</param>
/// <returns>运算成功返回true。如果向量ABC的元素个数不一致，运算失败，返回false。</returns>
bool vec_sub(vec A, vec B, vec C)
{
	real_t* A_data = A.data;
	real_t* B_data = B.data;
	real_t* C_data = C.data;
	if ((A.size == B.size) && (A.size == C.size))
	{
		for (int i = 0; i < A.size; i++) {
			C_data[i] = (A_data[i] - B_data[i]);
		}
		return true;
	}
	else {
		return false;
	}
}

/// <summary>
/// 实现向量与数的乘法
/// </summary>
/// <param name="A">向量A</param>
/// <param name="r">实数r</param>
/// <param name="B">向量B=A*r</param>
/// <returns>运算成功返回true。如果向量AB的元素个数不一致，运算失败，返回false。</returns>
bool vec_mul(vec A, real_t r, vec B)
{
	real_t* A_data = A.data;
	real_t* B_data = B.data;
	if (A.size == B.size)
	{
		for (int i = 0; i < A.size; i++) {
			B_data[i] = A_data[i] * r;
		}
		return true;
	}
	else {
		return false;
	}
}

/// <summary>
/// 计算两个向量之间的点乘。向量A和向量B的元素个数不需要一致，低维向量在高维度的分量默认为零。
/// </summary>
/// <param name="A">向量A</param>
/// <param name="B">向量B</param>
/// <returns>点乘的结果</returns>
real_t vec_dotProduct(vec A, vec B)
{
	real_t dotproduct = 0;
	real_t* A_data = A.data;
	real_t* B_data = B.data;
	int size = A.size;
	if (size < B.size)  size = B.size;

	for (int i = 0; i < size; i++) {
		dotproduct = dotproduct + (B_data[i] * A_data[i]);
	}
	return dotproduct;
}

/// <summary>
/// 计算向量模长。
/// </summary>
/// <param name="v">向量</param>
/// <returns>模长 = sqrtf(v[0]^2+v[1]^2+v[2]^2+…+v[n]^2) </returns>
real_t vec_length(vec v )
{
	real_t* data = v.data;
	real_t m = 0;
	for (int i = 0; i < v.size; i++) {
		m = m + (data[i] * data[i]);
	}
	return sqrtf(m);
}

/// <summary>
/// 计算向量模长的平方。
/// </summary>
/// <param name="v">向量</param>
/// <returns>模长平方 = (v[0]^2+v[1]^2+v[2]^2+…+v[n]^2) </returns>
real_t vec_length_square(vec v)
{
	real_t* data = v.data;
	real_t m = 0;
	for (int i = 0; i < v.size; i++) {
		m = m + (data[i] * data[i]);
	}
	return m;
}



/// <summary>
/// 计算矢量的曼哈顿距离Manhattan distance。
/// </summary>
/// <param name="v">向量</param>
/// <returns> 曼哈顿距离 = (|v[0]|+|v[1]|+|v[2]|+…+|v[n]|) </returns>
real_t vec_manhattan(vec v)
{
	real_t* data = v.data;
	real_t m = 0;
	for (int i = 0; i < v.size; i++) {
		m = m + ABS(data[i]);
	}
	return (m);
}


/// <summary>
/// 计算矢量v1到v2的曼哈顿距离Manhattan distance。
/// </summary>
/// <param name="v1">向量</param>
/// <param name="v2">向量</param>
/// <returns> 曼哈顿距离 = (|v1[0]-v2[0]|+|v1[1]-v2[1]|+……) </returns>
real_t vec_manhattan_2(vec v1,vec v2)
{
	real_t* data1 = v1.data;
	real_t* data2 = v2.data;
	real_t m = 0;
	real_t s;
	int size = v1.size > v2.size ? v2.size : v1.size;
	for (int i = 0; i < size; i++) {
		s = data1[i] - data2[i];
		m += ABS(s);
	}
	return (m);
}



/// <summary>
/// 随机初始化矩阵
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="seed">生成随机数的种子</param>
/// <param name="min">随机数的最小值</param>
/// <param name="max">随机数的最大值</param>
void mat_randomize(mat m, int seed, real_t min, real_t max)
{
	real_t** data = m.data;

	int meisen = 2147483647;   //m是梅森素数（形式为2^(p)-1的素数）：m=2^(31)-1  =  2147483647
	int a7_5 = 16807;		  //a=7^(5)
	unsigned long int x_i = seed;  //0-m之间的整数。是余数。
	real_t random_number;

	x_i = (a7_5 * x_i) % meisen;
	random_number = (real_t)x_i / meisen;  //现在random_number是0-1的分布	
	random_number = min + (max - min) * random_number; //把random_number缩放到min到max的分布

	for (int i = 0; i < m.rows; i++) {
		for (int j = 0; j < m.cols; j++) {
			x_i = (a7_5 * x_i) % meisen;
			random_number = (real_t)x_i / meisen;  //现在random_number是0-1的分布	
			random_number = min + (max - min) * random_number; //把random_number缩放到min到max的分布
			data[i][j] = random_number;
		}
	}
}


//实现矩阵的复制   A矩阵 => 矩阵B ，该矩阵复制给B的矩阵
bool mat_copyto(mat A,mat B)
{
	real_t** A_data = A.data;
	real_t** B_data = B.data;
	if ((A.cols == B.cols) && (A.rows == B.rows))
	{
		for (int i = 0; i < A.rows; i++) {
			for (int j = 0; j < A.cols; j++) {
				B_data[i][j] = A_data[i][j];
			}
		}
		return true;
	}
	else return false;
}

//实现矩阵的复制   A矩阵 <= 矩阵B ，B的矩阵复制给该矩阵
bool mat_copyfrom(mat A,mat B)
{
	real_t** A_data = A.data;
	real_t** B_data = B.data;
	if ((A.cols == B.cols) && (A.rows == B.rows))
	{
		for (int i = 0; i < A.rows; i++) {
			for (int j = 0; j < A.cols; j++) {
				A_data[i][j] = B_data[i][j];
			}
		}
		return true;
	}
	else return false;
}

/// <summary>
/// 初始化矩阵
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="value">矩阵元素初始化的值</param>
void mat_init(mat m, real_t value )
{
	real_t** data = m.data;
	int i, j;
	for (i = 0; i < m.rows; i++)
	{
		for (j = 0; j < m.cols; j++)
		{
			data[i][j] = value;
		}
	}
}

/// <summary>
/// 初始化区域。左上角为坐标原点，右下方是xy的正方向
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="left">区域左边（闭）</param>
/// <param name="top">区域上边（闭）</param>
/// <param name="right">区域右边（开）</param>
/// <param name="buttom">区域底边（开）</param>
/// <param name="value">初始化值</param>
void mat_initarea(mat m, int left, int top, int right, int buttom, real_t value)
{
	real_t** data = m.data;
	int i, j;
	if (left < 0)left = 0;
	if (top < 0)top = 0;
	if (right > m.cols)right = m.cols;
	if (buttom > m.rows)buttom = m.rows;
	for (i = top; i < buttom; i++)
	{
		for (j = left; j < right; j++)
		{
			data[i][j] = value;
		}
	}
}

/// <summary>
/// 对矩阵行矢量进行单位化
/// </summary>
/// <param name="m">矩阵</param>
void mat_row_vectorUnitization(mat m)
{
	int i, j;
	real_t** m_data = m.data;
	for (i = 0; i < m.rows; i++) {
		real_t d = 0;
		for (j = 0; j < m.cols; j++) {
			d = d + m_data[i][j] * m_data[i][j];
		}
		d = sqrtf(d);
		for (j = 0; j < m.cols; j++) {
			m_data[i][j] = m_data[i][j] / d;
		}
	}
}

/// <summary>
/// 对矩阵列矢量进行单位化
/// </summary>
/// <param name="m">矩阵</param>
void mat_col_vectorUnitization(mat m)
{
	int i, j;
	real_t** m_data = m.data;
	for (j = 0; j < m.cols; j++) {
		real_t d = 0;
		for (i = 0; i < m.rows; i++) {
			d = d + m_data[i][j] * m_data[i][j];
		}
		d = sqrtf(d);
		for (i = 0; i < m.rows; i++) {
			m_data[i][j] = m_data[i][j] / d;
		}
	}
}

/// <summary>
/// 单位化矩阵
/// </summary>
/// <param name="m">矩阵</param>
/// <returns>操作成功返回true。若不是方阵，操作失败，返回false</returns>
bool mat_eye(mat m)
{
	if (m.cols == m.rows)
	{
		mat_init(m, 0);
		for (int i = 0; i < m.rows; i++) {
			m.data[i][i] = 1;
		}
		return true;
	}
	else return false;
}

/// <summary>
/// 实现矩阵之间的加法
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="B">矩阵B</param>
/// <param name="C">矩阵C=A+B</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_add(mat A, mat B, mat C)
{
	real_t** A_data = A.data;
	real_t** B_data = B.data;
	real_t** C_data = C.data;
	if ((A.cols == B.cols) && (A.rows == B.rows) &&
		(A.cols == C.cols) && (A.rows == C.rows))
	{
		for (int i = 0; i < A.rows; i++) {
			for (int j = 0; j < A.cols; j++) {
				C_data[i][j] = (A_data[i][j] + B_data[i][j]);
			}
		}
		return true;
	}
	else {
		return false;
	}
}

/// <summary>
/// 实现矩阵之间的减法
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="B">矩阵B</param>
/// <param name="C">矩阵C=A-B</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_sub(mat A, mat B, mat C)
{
	real_t** A_data = A.data;
	real_t** B_data = B.data;
	real_t** C_data = C.data;
	if ((A.cols == B.cols) && (A.rows == B.rows) &&
		(A.cols == C.cols) && (A.rows == C.rows))
	{
		for (int i = 0; i < A.rows; i++) {
			for (int j = 0; j < A.cols; j++) {
				C_data[i][j] = (A_data[i][j] - B_data[i][j]);
			}
		}
		return true;
	}
	else {
		return false;
	}
}

/// <summary>
/// 实现矩阵之间的乘法
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="B">矩阵B</param>
/// <param name="C">矩阵C=A*B</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_mul(mat A, mat B, mat C)
{
	if ((A.cols == B.rows) && (A.rows == C.rows) && (B.cols == C.cols)) //先判断矩阵的尺寸能不能符合。
	{
		real_t** A_data = A.data;
		real_t** B_data = B.data;
		real_t** C_data = C.data;
		for (int i = 0; i < A.rows; i++) {
			for (int j = 0; j < B.cols; j++) {
				real_t value = 0;
				for (int k = 0; k < A.cols; k++) {
					value += (A_data[i][k] * B_data[k][j]);
				}
				C_data[i][j] = value;	//访问C.data[i][j]要先访问C,再移到data，再到data[i]，最后到data[i][j]，寻址效率比较低。所以用中间变量来提高速度。
			}
		}
		return true;
	}
	else {
		return false;
	}
}

/// <summary>
/// 实现矩阵元素之间的乘法
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="B">矩阵B</param>
/// <param name="C">矩阵C元素=A对应位置元素*B对应位置元素</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_data_mul(mat A, mat B, mat C)
{
	if ((A.cols == B.cols) && (A.rows == B.rows) && 
		(A.cols == C.cols) && (A.rows == C.rows)) //先判断矩阵的尺寸能不能符合。
	{
		real_t** A_data = A.data;
		real_t** B_data = B.data;
		real_t** C_data = C.data;
		for (int i = 0; i < A.rows; i++) {
			for (int j = 0; j < A.cols; j++) {
				C_data[i][j] = (A_data[i][j] * B_data[i][j]);
			}
		}
		return true;
	}
	else {
		return false;
	}
}

/// <summary>
/// 实现矩阵与向量的乘法
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="B">列矢B</param>
/// <param name="C">列矢C=矩阵A＊列矢B</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_mul_vec(mat A,vec B, vec C)
{
	if ((A.cols == B.size) && (A.rows == C.size)) //先判断矩阵的尺寸能不能符合。
	{
		real_t** A_data = A.data;
		real_t* B_data = B.data;
		real_t* C_data = C.data;

		for (int i = 0; i < A.rows; i++) {
			real_t value = 0;
			for (int j = 0; j < A.cols; j++) {
				value += (A_data[i][j] * B_data[j]);
			}
			C_data[i] = value;
		}
		return true;
	}
	else {
		return false;
	}
}

/// <summary>
/// 实现矩阵转置与向量的乘法
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="B">列矢B</param>
/// <param name="C">列矢C=(A^T)*B。等效式：行矢C=行矢B*矩阵A</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool matT_mul_vec(mat A,vec B, vec C)
{
	if ((A.rows == B.size) && (A.cols == C.size)) //先判断矩阵的尺寸能不能符合。
	{
		real_t** A_data = A.data;
		real_t* B_data = B.data;
		real_t* C_data = C.data;
		for (int i = 0; i < A.cols; i++) {
			real_t value = 0;
			for (int j = 0; j < A.rows; j++) {
				value += (A_data[j][i] * B_data[j]);
			}
			C_data[i] = value;	
		}
		return true;
	}
	else {
		return false;
	}
}

/// <summary>
/// 实现矩阵与数的乘法
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="r">实数r</param>
/// <param name="B">矩阵B</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_mul_real(mat A, real_t r, mat B)
{
	if ((A.cols == B.cols) && (A.rows == B.rows))
	{
		real_t** A_data = A.data;
		real_t** B_data = B.data;
		for (int i = 0; i < A.rows; i++) {
			for (int j = 0; j < A.cols; j++) {
				B_data[i][j] = A_data[i][j] * r;
			}
		}
		return true;
	}
	else {
		return false;
	}
}

/// <summary>
/// 矢量乘法生成矩阵。矩阵A=列矢B＊列矢C^T
/// </summary>
/// <param name="A">矩阵A=列矢B＊列矢C^T</param>
/// <param name="B">列矢B</param>
/// <param name="C">列矢C</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_vec_mul_vecT(mat A, vec B, vec C)
{
	if ((A.rows == B.size) && (A.cols == C.size)) //先判断矩阵的尺寸能不能符合。
	{
		real_t** A_data = A.data;
		real_t* B_data = B.data;
		real_t* C_data = C.data;

		for (int i = 0; i < A.rows; i++) {
			for (int j = 0; j < A.cols; j++) {
				A_data[i][j] = B_data[i] * C_data[j];
			}
		}
		return true;
	}
	else {
		return false;
	}
}



///////////////////////////////////////////////////////////////////

/// <summary>
/// 高斯消元法求逆矩阵。
/// 如果运算成功，原来的矩阵A变成单位矩阵。
/// 如果运算失败，原来的矩阵A也会被破坏。最好事先拷贝一份副本。
/// </summary>
/// <param name="A">原来的矩阵，成功求逆后变成单位矩阵</param>
/// <param name="iA">逆矩阵，运算结束后变成原来矩阵A的逆</param>
/// <param name="precision">运算精度</param>
/// <returns>如果是满秩的方阵，可以求逆，返回true。否则无法正确求逆，返回false</returns>
bool  mat_inversion(mat A, mat iA, real_t precision )
{
	//方阵满秩才能求逆
	if ((A.rows == A.cols) && (A.rows == iA.rows) && (A.cols == iA.cols)) {  //检查行数和列数是否匹配
		real_t** A_data = A.data;
		real_t** iA_data = iA.data;
		mat_eye(iA);   //先把iA初始化。
		int i = 0, j = 0, k = 0;
		real_t d;
		real_t* line_pointer;		//行指针，用于交换行向量时使用。

		//用列主元高斯消元把A化成上三角形式。
		for (j = 0; j < A.cols - 1; j++) {		//从第一列开始，逐列消去下三角区元素。
			d = 0;
			k = j;
			//选择列主元
			for (i = j; i < A.rows; i++) {
				if (ABS(A_data[i][j]) > d) {
					d = ABS(A_data[i][j]);
					k = i;		//记录列主元行数。
				}
			}
			if (d < precision) return false;
			//如果主元等于0（绝对值小于precision），说明矩阵值为零，不能求逆。这时可以通过QR分解求伪逆。
			if (k != j) {
				//如果k!=j, 说明需要交换行以选择最大主元。只需要交换行指针即可，可节省计算量。
				line_pointer = A_data[j];
				A_data[j] = A_data[k];
				A_data[k] = line_pointer;
				line_pointer = iA_data[j];
				iA_data[j] = iA_data[k];
				iA_data[k] = line_pointer;
			}
			//高斯消元。
			for (i = j + 1; i < A.rows; i++) {
				d = A_data[i][j] / A_data[j][j];		//计算行向量相减的倍数
				for (k = j; k < A.cols; k++) {
					A_data[i][k] = A_data[i][k] - d * A_data[j][k];
					iA_data[i][k] = iA_data[i][k] - d * iA_data[j][k];		//iA也要跟着变换。
				}
				for (k = 0; k < j; k++) {
					iA_data[i][k] = iA_data[i][k] - d * iA_data[j][k];
				}
			}
		}
		if (ABS(A_data[A.rows - 1][A.rows - 1]) < precision) {
			return false;			//前面检查最大主元绝对值大于零时，没有检查最后一个主元。
		}
		//现在A是上三角形式，进一步把A上三角区域消去。
		for (j = A.cols - 1; j >= 0; j--) {			//从最右边的一列开始往左消
			for (i = j - 1; i >= 0; i--) {
				d = A_data[i][j] / A_data[j][j];			//计算行向量相减的倍数
				A_data[i][j] = 0;
				for (k = 0; k < A.cols; k++) {
					iA_data[i][k] = iA_data[i][k] - d * iA_data[j][k];  //变换iA。
				}
			}
		}

		//现在A是对角的，但不是单位的。下面把A化成单位的。
		for ( i = 0; i < A.rows; i++) {
			for ( j = 0; j < A.cols; j++) {
				iA_data[i][j] = iA_data[i][j] / A_data[i][i];
			}
			A_data[i][i] = 1;
		}

		return true;
	}
	else
		return false;
}

/// <summary>
/// 计算方阵的值。按照高斯消元法将矩阵分解成上三角矩阵，然后计算对角线积。
/// 如果运算成功，原来的矩阵m变成上三角矩阵。
/// 如果运算失败，原来的矩阵m会被破坏。最好事先拷贝一份副本。
/// </summary>
/// <param name="m">原来的矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>矩阵的值。如果矩阵是满秩的方阵，则可以计算值。否则返回0.</returns>
real_t mat_det(mat m,real_t precision )
{
	if (m.rows == m.cols) {  //检查行数和列数是否匹配
		// Matrix A(m.rows, m.cols, m.buffer);    	
		real_t** A_data = m.data;   //

		int i = 0, j = 0, k = 0;
		int sign = 1;   //记录行交换的次数。以计算值的正负号。
		real_t d;
		real_t* line_pointer;		//行指针，用于交换行向量时使用。

		//用列主元高斯消元把A化成上三角形式。
		for (j = 0; j < m.cols - 1; j++) {		//从第一列开始，逐列消去下三角区元素。
			d = 0;
			k = j;
			//选择列主元
			for (i = j; i < m.rows; i++) {
				if (ABS(A_data[i][j]) > d) {
					d = ABS(A_data[i][j]);
					k = i;		//记录列主元行数。
				}
			}
			if (ABS(d) < precision) return 0;   //如果主元等于0（绝对值小于precision），说明矩阵值为0。
			if (k != j) {	//如果k!=j, 说明需要交换行以选择最大主元。只需要交换行指针即可，可节省计算量。
				line_pointer = A_data[j];
				A_data[j] = A_data[k];
				A_data[k] = line_pointer;
				sign = -sign;  //每次行交换，都会改变一次正负号。
			}
			//高斯消元。
			for (i = j + 1; i < m.rows; i++) {
				d = A_data[i][j] / A_data[j][j];		//计算行向量相减的倍数
				for (k = j; k < m.cols; k++) {
					A_data[i][k] = A_data[i][k] - d * A_data[j][k];
				}
			}
		}
		if (ABS(A_data[m.rows - 1][m.rows - 1]) < precision) {
			return 0;			//前面检查最大主元绝对值大于零时，没有检查最后一个主元。
		}

		//现在A是上三角形式，其对角元素的乘积即为矩阵的值。
		d = A_data[0][0];
		for (i = 1; i < m.rows; i++) {			//从最右边的一列开始往左消
			d = d * A_data[i][i];
		}
		d = d * sign;
		return d;
	}
	else
		return 0;
}

/// <summary>
/// 计算矩阵的范数，即矩阵每个元素的平方和
/// </summary>
/// <param name="m">矩阵</param>
/// <returns>矩阵所有元素的平方和</returns>
real_t mat_norm(mat m) {
	real_t r = 0;
	real_t** dat = m.data;
	for (int i = 0; i < m.rows; i++) {
		for (int j = 0; j < m.cols; j++) {
			r += dat[i][j] * dat[i][j];
		}
	}
	return r;
}

/// <summary>
/// 实现矩阵转置
/// </summary>
/// <param name="m">矩阵</param>
/// <returns>如果矩阵是方阵，返回true，矩阵m变成转置矩阵。否则返回false，矩阵m不变</returns>
bool mat_transpose(mat m)
{
	if (m.cols == m.rows)
	{
		real_t temp; //临时变量
		for (int i = 0; i < m.rows; i++) {
			for (int j = i; j < m.cols; j++) {  
				//注意，这里的循环要从j=i开始，否则会执行两次转置操作，相互抵消。
				temp = m.data[i][j];
				m.data[i][j] = m.data[j][i];
				m.data[j][i] = temp;
			}
		}
		return true;
	}
	else {
		return false;
	}
}

/// <summary>
/// 实现矩阵转置
/// </summary>
/// <param name="A">原矩阵</param>
/// <param name="AT">转置后的矩阵</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_transpose_to(mat A, mat AT)
{
	if ((AT.cols == A.rows) && (AT.rows == A.cols))
	{
		for (int i = 0; i < A.rows; i++) {
			for (int j = 0; j < A.cols; j++) {
				AT.data[j][i] = A.data[i][j];
			}
		}
		return true;
	}
	else return false;
}
 
/// <summary>
/// 判断矩阵是否相等
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="B">矩阵B</param>
/// <param name="precision">运算精度</param>
/// <returns>A等于B，返回true。否则返回false</returns>
bool mat_isEqual(mat A, mat B, real_t precision )
{
	if ((B.rows == A.rows) && (B.cols == A.cols))
	{
		for (int i = 0; i < A.rows; i++) {
			for (int j = 0; j < A.cols; j++) {
				if ( (B.data[i][j] - A.data[i][j])* (B.data[i][j] - A.data[i][j]) > precision) return false;
			}
		}
		return true;
	}
	else {
		return false;
	}
}

/// <summary>
/// 判断该矩阵是否是上三角矩阵。只看下三角区域是否全为0，不管是不是方针。
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>m是上三角矩阵返回true，否则返回false</returns>
bool mat_isUpperTriangle(mat m,real_t precision )
{
	{
		real_t** m_data = m.data;
		for (int i = 0; i < m.rows; i++) {
			for (int j = 0; j < i && j < m.cols; j++) {
				if (m_data[i][j] > precision) return false;
			}
		}
		return true;
	}
	// return false;
}

/// <summary>
/// 判断该矩阵是否是下三角矩阵。只看上三角区域是否全为0，不管是不是方针。
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>m是下三角矩阵返回true，否则返回false</returns>
bool mat_isLowerTriangle(mat m,real_t precision )
{
//	if (m.cols == m.rows)
	{
		real_t** m_data = m.data;
		for (int i = 0; i < m.rows; i++) {
			for (int j = i + 1; j < m.cols; j++) {
				if (m_data[i][j] > precision) return false;
			}
		}
		return true;
	}
//	return false;
}

/// <summary>
/// 判断该矩阵是否是对称矩阵
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>m是对称矩阵返回true，否则返回false</returns>
bool mat_isSymmetric(mat m,real_t precision )
{
	if (m.cols == m.rows)
	{
		real_t** m_data = m.data;
		for (int i = 0; i < m.rows; i++) {
			for (int j = i + 1; j < m.cols; j++) {
				if (ABS(m_data[j][i] - m_data[i][j]) > precision) return false;
			}
		}
		return true;
	}
	else return false;
}

/// <summary>
/// 判断该矩阵是否单位矩阵
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>m是单位矩阵返回true，否则返回false</returns>
bool mat_isEye(mat m,real_t precision )
{
	if (m.cols == m.rows)
	{
		real_t** m_data = m.data;
		//检查非对角元素是否为零
		for (int i = 0; i < m.rows; i++) {
			for (int j = i + 1; j < m.cols; j++) {
				if ((m_data[i][j]) * (m_data[i][j]) > precision) return false;
				if ((m_data[j][i]) * (m_data[j][i]) > precision) return false;
			}
		}
		//检查对角元素是否为1.
		for (int i = 0; i < m.rows; i++) {
			if ((m_data[i][i] - 1) * (m_data[i][i] - 1) > precision) return false;
		}
		return true;
	}
	else return false;
}

/// <summary>
/// 判断该矩阵 对角线是否0或1，其他元素都为0
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>m是单位矩阵返回true，否则返回false</returns>
bool mat_isDiagonal_01(mat m, real_t precision)
{
	if (m.cols == m.rows)
	{
		real_t** m_data = m.data;
		//检查非对角元素是否为零
		for (int i = 0; i < m.rows; i++) {
			for (int j = i + 1; j < m.cols; j++) {
				if (ABS(m_data[i][j]) > precision) return false;
				if (ABS(m_data[j][i]) > precision) return false;
			}
		}
		//检查对角元素是否为1.
		for (int i = 0; i < m.rows; i++) {
			if (ABS(m_data[i][i] - 1) < precision || ABS(m_data[i][i]) < precision) {

			}
			else			
				return false;

		}
		return true;
	}
	else return false;
}




///////////////////  下面算是中级内容了     ////////////////////////////



/// <summary>
/// 豪斯霍尔德变换。
/// 变换后，A = H * A 。H是豪斯霍尔德矩阵，以豪斯霍尔德矢量ｖ和参数beta的形式出现。 H=I-(1/beta)v*v^(T)
/// 变换后的A的x1列，y1行以下的元素是0元素，把x1列y1行开始的列向量变成(a,0,0,...)^(T)形式的向量。得到的豪斯霍尔德矩阵以向量形式保存在v中。
/// 注意事项：注意数组v的大小，防止越界。
/// 小于x1的列是不变换的。
/// </summary>
/// <param name="A">矩阵</param>
/// <param name="x1">矩阵豪斯霍尔德变换的起始列</param>
/// <param name="y1">矩阵豪斯霍尔德变换的起始行</param>
/// <param name="v">豪斯霍尔德矢量</param>
/// <param name="beta">豪斯霍尔德矢量的参数</param>
/// <param name="isH_ensured">是否已经确认豪斯霍尔德矢量，如果未确认，函数会先计算更新豪斯霍尔德矢量和参数</param>
/// <param name="precision">计算精度</param>
void mat_Householder_left(mat A, int x1, int y1, real_t* v, real_t* beta, bool isH_ensured , real_t precision )
{
	int i, j;
	real_t d;
	real_t** A_data = A.data;
	if (isH_ensured == false) {  //如果还没确定Householder向量，先计算Householder向量。
		d = 0;
		//计算矩阵A第x1列y1行开始的列向量模长
		for (i = y1; i < A.rows; i++) {
			d += A_data[i][x1] * A_data[i][x1];
		}
		d = sqrtf(d);   //现在d是模长。
		//计算向量v。向量v有A.rows-1-y1个元素。
		(beta[0]) = d * (d - A_data[y1][x1]);
		//如果beta等于零，则H为单位矩阵，向量v为零向量。(从beta的计算可知，beta是大于零的)
		if (beta[0] < precision) {
			beta[0] = 1;   	//β不为零是防止下面除法运算出错。 
			for (i = y1; i < A.rows; i++) {
				v[i - y1] = 0;
			}
		}
		else {				//否则更新Hi。同时更新矩阵A的x1列
			v[0] = A_data[y1][x1] - d;
			A_data[y1][x1] = d;
			for (i = y1 + 1; i < A.rows; i++) {
				v[i - y1] = A_data[i][x1];
				A_data[i][x1] = 0;
			}
		}

		real_t div_beta = 1 / (beta[0]);
		//计算A左乘Householder矩阵 。A=H(i)*A。从x1+1列开始变换。
		for (j = x1 + 1; j < A.cols; j++) {	// 
			d = 0;
			for (i = y1; i < A.rows; i++) {	//计算Householder向量和A列矢的点积
				d = d + v[i - y1] * A_data[i][j];
			}
			d = d * div_beta;
			for (i = y1; i < A.rows; i++) {
				A_data[i][j] = A_data[i][j] - d * v[i - y1];
			}
		}
	}
	else {
		real_t div_beta = 1 / (beta[0]);
		//计算A左乘Householder矩阵 。A=H(i)*A。从x1列开始变换。
		for (j = x1; j < A.cols; j++) {
			d = 0;
			for (i = y1; i < A.rows; i++) {	//计算Householder向量和A列矢的点积
				d = d + v[i - y1] * A_data[i][j];
			}
			d = d * div_beta;
			for (i = y1; i < A.rows; i++) {
				A_data[i][j] = A_data[i][j] - d * v[i - y1];
			}
		}
	}
}

/// <summary>
/// 豪斯霍尔德变换。
/// 变换后，A = A * H。H是豪斯霍尔德矩阵，以豪斯霍尔德矢量ｖ和参数beta的形式出现。 H=I-(1/beta)v*v^(T)
/// 变换后的A的x1列，y1行以右的元素是0元素，把x1列y1行开始的行向量变成(a,0,0,...)形式的向量。得到的豪斯霍尔德矩阵以向量形式保存在v中。
/// 注意事项：注意数组v的大小，防止越界。
/// 小于y1的行是不变换的。
/// </summary>
/// <param name="A">矩阵</param>
/// <param name="x1">矩阵豪斯霍尔德变换的起始列</param>
/// <param name="y1">矩阵豪斯霍尔德变换的起始行</param>
/// <param name="v">豪斯霍尔德矢量</param>
/// <param name="beta">豪斯霍尔德矢量的参数</param>
/// <param name="isH_ensured">是否已经确认豪斯霍尔德矢量，如果未确认，函数会先计算更新豪斯霍尔德矢量和参数</param>
/// <param name="precision">计算精度</param>
void mat_Householder_right(mat A, int x1, int y1, real_t* v, real_t* beta, bool isH_ensured , real_t precision )
{
	int i, j;
	real_t d;
	real_t** A_data = A.data;
	if (isH_ensured == false) {
		//先求Householder矩阵的向量
		d = 0;
		//计算矩阵A第x1列y1行开始的行向量模长
		for (i = x1; i < A.cols; i++) {
			d += A_data[y1][i] * A_data[y1][i];
		}
		d = sqrtf(d);   //现在d是模长。
		//计算向量v。向量v有A.rows-1-y1个元素。
		beta[0] = d * (d - A_data[y1][x1]);
		//如果beta等于零，则H为单位矩阵，向量v为零向量。(从beta的计算可知，beta是大于零的)
		if (beta[0] < precision) {
			beta[0] = 1;   	//β不为零是防止下面除法运算出错。 
			for (i = x1; i < A.cols; i++) {
				v[i - x1] = 0;
			}
		}
		else {				//否则更新Hi。同时更新矩阵A的y1行
			v[0] = A_data[y1][x1] - d;
			A_data[y1][x1] = d;
			for (i = x1 + 1; i < A.cols; i++) {
				v[i - x1] = A_data[y1][i];
				A_data[y1][i] = 0;
			}
		}

		real_t div_beta = 1 / (beta[0]);
		//计算A右乘Householder矩阵 。A=A*H(i)。从y1+1行开始变换。
		for (j = y1 + 1; j < A.rows; j++) {
			d = 0;
			for (i = x1; i < A.cols; i++) {			//计算Householder向量和A列矢的点积
				d = d + v[i - x1] * A_data[j][i];
			}
			d = d * div_beta;
			for (i = x1; i < A.cols; i++) {
				A_data[j][i] = A_data[j][i] - d * v[i - x1];
			}
		}
	}
	else {
		real_t div_beta = 1 / (beta[0]);
		//计算A右乘Householder矩阵 。A=A*H(i)。从y1行开始变换。
		for (j = y1; j < A.rows; j++) {
			d = 0;
			for (i = x1; i < A.cols; i++) {			//计算Householder向量和A列矢的点积
				d = d + v[i - x1] * A_data[j][i];
			}
			d = d * div_beta;
			for (i = x1; i < A.cols; i++) {
				A_data[j][i] = A_data[j][i] - d * v[i - x1];
			}
		}
	}
}



/// <summary>
/// 使用子函数的豪斯霍尔德QR分解。
/// 不要求m矩阵是方阵。A尺寸与m矩阵一样。H是和行数一样的方阵。
/// 不管A是不是满秩，分解后的H是满秩的。
/// </summary>
/// <param name="m">原矩阵=H*A</param>
/// <param name="operation_buffer">运算的缓存矩阵。m.cols行，m.rows+1列。要cols个豪斯霍尔德向量，第一个向量rows个元素。每个向量还要一个beta参数</param>
/// <param name="H">分解的正交矩阵</param>
/// <param name="A">分解的上三角矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_QR_Householder(mat m,  mat H, mat A, real_t precision ) //m矩阵=H*A.
{
	if (m.cols < 1)return true;  //如果列数小于1，就没必要运算了。

	if ((m.rows == H.cols) && (m.rows == H.rows) &&   
		(m.cols == A.cols) && (m.rows == A.rows))
	{
		int  j;
		int n = m.cols >= m.rows ? (m.rows - 1) : (m.cols - 1);//选行列数比较小的值。

		//要cols个豪斯霍尔德向量，第一个向量rows个元素，第二个向量rows-1个元素.......
		//第i个向量rows-(i-1)个元素，最后一个向量rows-(cols-1-1)个元素。
		// 每个豪斯霍尔德向量需要一个beta参数。
		mat operation_buffer;
		vec beta_buffer;
		real_t** v;
		real_t beta[1];
		mat_allocate(&operation_buffer, m.cols, m.rows);
		vec_allocate(&beta_buffer, m.cols);
		v = operation_buffer.data;
		//内存分配完成。	
		mat_copyto(m, A);//把该矩阵复制给A

		//下面开始变换
//		mat_print(A);
		for (j = 0; j <= n; j++) {
			mat_Householder_left(A, j, j, v[j], beta, false, precision);
			beta_buffer.data[j] = beta[0];
			//mat_print(A);
		}
		//printfs("\r\n  打印矩阵mat_print(A);完成  +++++++++++++++++++ \r\n";

		//下面计算H，H=H(0)*H(1)*……*H(n-2)*I。

		//初始化H为单位矩阵
		mat_eye(H);
		//开始计算H
		for (j = n; j >= 0; j--) {
			beta[0] = beta_buffer.data[j];
			mat_Householder_left(H, j, j, v[j], beta, true, precision);
		}
		//mat_print(H);
		//mat_print(operation_buffer);
		//mat_print(A);

		mat_deallocate(&operation_buffer);
		vec_deallocate(&beta_buffer);

		return true;
	}
	else return false;
}


/// <summary>
/// 利用Givens变换实现QR分解。当该矩阵下三角区域存在大量零元的时候，Givens变换进行QR分解速度更快。
///  G是m.rows大小的方针，A的矩阵尺寸和m的矩阵一致。
/// </summary>
/// <param name="m">原矩阵=G*A</param>
/// <param name="G">分解的正交矩阵</param>
/// <param name="A">分解的上三角矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_QR_Givens(mat m, mat G, mat A, real_t precision )  //m = G*A。G是正交矩阵，A是上三角矩阵。
{
	if (	// G是m.rows大小的方针，  m的矩阵和A的矩阵尺寸一致。
		(m.rows == G.cols) && (m.rows == G.rows) &&
		(m.cols == A.cols) && (m.rows == A.rows))
	{
		int i, j, i2, j2;
		real_t sin, cos; //保存吉文斯矩阵参数
		real_t d, d1, d2;	 //d1,d2用来保存矩阵在进行吉文斯变换时的中间参数
		real_t** A_data = A.data;
		real_t** G_data = G.data;
		//先初始化G为单位矩阵，和复制该矩阵到A
		mat_eye(G);
		mat_copyto(m, A);
		//开始Givens变换
		{
			for (j = 0; j < m.cols; j++) {		//逐列消元
				for (i = j + 1; i < m.rows; i++) {  //消下三角区域元素
					if (ABS(A_data[i][j]) >= precision) {  //只有下三角元素不为零才需要吉文斯变换

						//计算吉文斯矩阵参数
						d = sqrtf(A_data[j][j] * A_data[j][j] + A_data[i][j] * A_data[i][j]);
						cos = A_data[j][j] / d;    //G= |cos,  sin|  ,G是单位正交对称矩阵，G^(-1)=G
						sin = A_data[i][j] / d;	//		|sin, -cos|

						//更新A：A左乘G(i).  A=G(i)*A
						for (j2 = j; j2 < m.cols; j2++) {    //从第j列开始计算。因为0--(j-1)列下三角都是0了，不用去计算。
							d1 = A_data[j][j2];
							d2 = A_data[i][j2];
							A_data[j][j2] = d1 * cos + d2 * sin;
							A_data[i][j2] = d1 * sin - d2 * cos;
						}
						//更新G：G右乘G(i).  G=G*G(i)  
						for (i2 = 0; i2 < i + 1; i2++) {		//有优化空间。可以在草稿纸上记录零元区域的变换情况，从而避免不必要的零元计算。
							d1 = G_data[i2][j];						//零元区域变化是有规律的。这里已经避免了零元区域的计算了。
							d2 = G_data[i2][i];
							G_data[i2][j] = d1 * cos + d2 * sin;
							G_data[i2][i] = d1 * sin - d2 * cos;
						}
					}
				}
			}
		}
		return true;
	}
	else return false;
}

/// <summary>
/// 使用改进“格拉姆施密特过程”的算法实现QR分解。
/// 格拉姆施密特过程不要求原矩阵是方阵，只要Q的行列数和原矩阵m一致，R是与列数相同的方阵。
/// 只有对角元全为正的时候，QR分解才是唯一的。不同方法的QR分解，分解出的Q，R不同行列向量正负号可能不同。Q的秩与原矩阵相同。
/// </summary>
/// <param name="m">原矩阵=Q*R</param>
/// <param name="Q">分解的正交矩阵</param>
/// <param name="R">分解的上三角矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_QR_Schmidt(mat m, mat Q, mat R, real_t precision )  //把该矩阵分解成Q和R（单位正交矩阵和上三角矩阵）
{
	if ((m.cols == Q.cols) && (m.rows == Q.rows) &&    //格拉姆施密特过程不要求Q是方阵，只要求Q与原矩阵尺寸一样
		(m.cols == R.cols) && (m.cols == R.rows))	   //R是与原矩阵列数相同的上三角方阵。
	{
		real_t** Q_data = Q.data;
		real_t** R_data = R.data;

		int i, j;//循环行和列的中间变量

		//把该矩阵复制给Q，把R初始化为单位矩阵
		mat_copyto(m, Q);
		mat_eye(R);
		//开始施密特正交化
		{
			real_t r;
			for (j = 0; j < Q.cols; j++)
			{
				//计算列向量模长
				r = 0;
				for (i = 0; i < Q.rows; i++)
				{
					r += Q_data[i][j] * Q_data[i][j];
				}
				if (r <= precision ) // 如果模长小于精确度,就认为它等于0。模长本来就大于0，不必取绝对值
				{
					for (i = 0; i < Q.rows; i++)
					{
						Q_data[i][j] = 0; //如果该矢量模长为零，则所有元素为零。
					}
					R_data[j][j] = 1;		//R对应的该列设为1。R可以取任何值，为了让R满秩，对角元素设为1.
				}
				else  //如果模长不为零，进行单位化
				{
					r = sqrtf(r);
					for (i = 0; i < Q.rows; i++)
					{
						Q_data[i][j] = Q_data[i][j] / r; //单位化该列向量
					}
					R_data[j][j] = r;
				}
				int j2; //下面更新矩阵Q和R，把Q的其他列向量减去q1的分量。要用到‘第二个循环中’列数的中间变量
				for (j2 = j + 1; j2 < Q.cols; j2++)
				{
					r = 0;
					for (i = 0; i < Q.rows; i++)
					{
						r += Q_data[i][j] * Q_data[i][j2];  //计算其他向量投影在Q_data(i)[j]上的长度。注：Q_data(i)[j]表示j列的向量
					}
					for (i = 0; i < Q.rows; i++)
					{
						Q_data[i][j2] = Q_data[i][j2] - (r * Q_data[i][j]);  //更新矩阵Q，使其他列向量消去Q_data(i)[j]的分量。
					}
					R_data[j][j2] = r;	//更新矩阵R
				}
			}
		}
		return true;
	}
	else return false;
}

/// <summary>
/// 线性方程组最小二乘法求解。m*X=mX，已知m和mX，求最佳的X。
/// 这里使用豪斯霍尔德变换进行QR分解，计算最小二乘。如果矩阵中间有零元行列，会算错。见下面代码注释。
/// </summary>
/// <param name="m">矩阵m</param>
/// <param name="operation_buffer">运算的缓存矩阵。m.rows*2+m.cols+1行，max(m.rows+1,m.cols)列。</param>
/// <param name="X">待求解的矢量X</param>
/// <param name="mX">已知的矢量mX</param>
/// <param name="is_resetX">是否重置矢量X。如果重置,计算X前会初始化为零矢量</param>
/// <param name="precision">运算精度</param>
/// <returns>操作是否成功。若矩阵不匹配，或者矩阵不是满秩，运算失败，返回false</returns>
bool mat_LinearEquation_LinearLeastSquare_Househoulder(mat m, vec X,vec mX,bool is_resetX, real_t precision )
{
	if ((m.cols == X.size) && (m.rows == mX.size)) {  //先验证矩阵m和矢量X，mX的运算是否合法。

		//X有公式解： X = [invers(mt*m)]*mT*mX.      invers表示求逆。mt表示m的转置。
		//QR分解： 令m=q*r,mt=rt*qt，因为qt*q=I，所以X=[invers(mt*m)]*mT*mX=[invers(rt*r)]*rt*qt*mX。

		// (m.rows, m.cols);
		mat R; 
		// (m.rows, m.rows);
		mat Q;	
		// (1,m.rows)
		vec QtB;

		mat_allocate(&R, m.rows, m.cols);
		mat_allocate(&Q, m.rows, m.rows);
		vec_allocate(&QtB, m.rows);
		 
		real_t** R_data = R.data;
		real_t* QtB_data = QtB.data;
		real_t* X_data = X.data;

		real_t d;
		int i, j;
		mat_QR_Householder(m, Q, R, precision);		//m=Q*R,==> R*X=Q^(T)*mX。

		matT_mul_vec(Q, mX, QtB);

		//mat_print(m);
		//mat_print(Q);
		//mat_print(R);
		//vec_print(QtB);

		if (is_resetX) {
			vec_init(X, 0);  //先初始化待求矢量为零矢量。后面计算的时候，X中没计算的元素就默认为零元素。
		}

		//下面用R和QtB回代X.   //R*X=QtB。 //R可能不是下三角。
		if (R.rows >= R.cols) {  //行数大于列数。

			//注意，行数大于列数的情况比较特别。举个例子
			// R= | r1,r2 | X= (x1,x2), Qtb= ( q1,q2,q3,q4 )
			//    | 0, r3 |
			//    | 0, 0  |
			//    | 0, 0  |
			//这里的Qtb的q3,q4如果不考虑，计算的结果有问题。
			
		//	mat_print(R);

			for (i = R.cols - 1; i >= 0; i--) {
				d = 0;
				for (j = i + 1; j < R.cols; j++) {
					d = d + R_data[i][j] * X_data[j];
				}
				if (ABS(R_data[i][i]) > precision) {
					//检查R的对角元素是否零元。如果不是，可以直接求x。否则x不唯一，则令该行x元素为0 。
					X_data[i] = (QtB_data[i] - d) / R_data[i][i];
				}
				else {
					X_data[i] = 0;
				}
			}

			//注意，行数大于列数的情况比较特别。下面的例子更特殊：
			// R= | r1,0,r2,r4  |      X= (x1,x2,x3,x4), Qtb= ( q1,q2,q3,q4 )
			//    | 0, 0,r3,r5  |
			//    | 0, 0 ,r3,r6  |
			//    | 0, 0 ,0 ,r7  |
			//这里矩阵R的子矩阵：
			//    | r3,r5  |  对应的 X的参数,（x2,x3,x4），  （q2,q3,q4 ）还要再用一次QR分解，才能解出X。
			//    | r3,r6  |
			//    | 0 ,r7  |
			//所以Household的QR分解最小二乘，有缺陷。
		}
		else {
			//注意，列数大于列行的情况。举个例子
			// R= | r1,r2,r3 | X= (x1,x2,x3), Qtb= ( q1,q2 )
			//    | 0, r4,r5 |
			//这里先从x2算起。x3用默认的值。

			for (i = R.rows - 1; i >= 0; i--) {
				d = 0;
				for (j = i + 1; j < R.cols; j++) {
					d = d + R_data[i][j] * X_data[j];
				}
				if (ABS(R_data[i][i]) > precision) {
					//检查R的对角元素是否零元。如果不是，可以直接求x。否则x不唯一，则令该行x元素为0 。
					X_data[i] = (QtB_data[i] - d) / R_data[i][i];
				}
				else {
					X_data[i] = 0;
				}
			}
		}

		//mat_print(Q);
		//mat_print(R);
		//mat_mul(Q, R, m);
		//mat_print(m);

		//vec_print(X);
		//vec_print(mX);

		mat_deallocate(&Q);
		mat_deallocate(&R);
		vec_deallocate(&QtB);

		return true;
	}
	else return false;
}

/// <summary>
/// 线性方程组最小二乘法求解。m*X=mX，已知m和mX，求最佳的X。
/// 这里使用Schmidt变换进行QR分解，计算最小二乘。 这个算法更普适。
/// </summary>
/// <param name="m">矩阵m</param>
/// <param name="operation_buffer">运算的缓存矩阵。m.rows*2+m.cols+1行，max(m.rows+1,m.cols)列。</param>
/// <param name="X">待求解的矢量X</param>
/// <param name="mX">已知的矢量mX</param>
/// <param name="is_resetX">是否重置矢量X。如果重置,计算X前会初始化为零矢量</param>
/// <param name="precision">运算精度</param>
/// <returns>操作是否成功。若矩阵不匹配，或者矩阵不是满秩，运算失败，返回false</returns>
bool mat_LinearEquation_LinearLeastSquare(mat m, vec X, vec mX, bool is_resetX, real_t precision)
{
	if ((m.cols == X.size) && (m.rows == mX.size)) {  //先验证矩阵m和矢量X，mX的运算是否合法。

		//X有公式解： X = [invers(mt*m)]*mT*mX.      invers表示求逆。mt表示m的转置。
		//QR分解： 令m=q*r,mt=rt*qt，因为qt*q=I，所以X=[invers(mt*m)]*mT*mX=[invers(rt*r)]*rt*qt*mX。

		// (m.rows, m.cols);
		mat R;
		// (m.rows, m.cols);
		mat Q;
		// (1,m.rows)
		vec QtB;

		mat_allocate(&R, m.cols, m.cols);
		mat_allocate(&Q, m.rows, m.cols);
		vec_allocate(&QtB, m.cols);

		real_t** R_data = R.data;
		real_t* QtB_data = QtB.data;
		real_t* X_data = X.data;

		real_t d;
		int i, j;
		mat_QR_Schmidt(m, Q, R, precision);		//m=Q*R,==> R*X=Q^(T)*mX。

		matT_mul_vec(Q, mX, QtB);

		//mat_print(m);
		//mat_print(Q);
		//mat_print(R);
		//vec_print(QtB);

		if (is_resetX) {
			vec_init(X, 0);  //先初始化待求矢量为零矢量。后面计算的时候，X中没计算的元素就默认为零元素。
		}

		//下面用R和QtB回代X.   //R*X=QtB。 //R可能不是下三角。
		if (R.rows >= R.cols) {  //行数大于列数。

			//注意，行数大于列数的情况比较特别。举个例子
			// R= | r1,r2 | X= (x1,x2), Qtb= ( q1,q2,q3,q4 )
			//    | 0, r3 |
			//    | 0, 0  |
			//    | 0, 0  |
			//这里的Qtb的q3,q4如果不考虑，计算的结果有问题。

		//	mat_print(R);

			for (i = R.cols - 1; i >= 0; i--) {
				d = 0;
				for (j = i + 1; j < R.cols; j++) {
					d = d + R_data[i][j] * X_data[j];
				}
				if (ABS(R_data[i][i]) > precision) {
					//检查R的对角元素是否零元。如果不是，可以直接求x。否则x不唯一，则令该行x元素为0 。
					X_data[i] = (QtB_data[i] - d) / R_data[i][i];
				}
				else {
					X_data[i] = 0;
				}
			}

		}
		else { //施密特正交分解不会出现这个情况。
			//注意，列数大于列行的情况。举个例子
			// R= | r1,r2,r3 | X= (x1,x2,x3), Qtb= ( q1,q2 )
			//    | 0, r4,r5 |
			//这里先从x2算起。x3用默认的值。

			for (i = R.rows - 1; i >= 0; i--) {
				d = 0;
				for (j = i + 1; j < R.cols; j++) {
					d = d + R_data[i][j] * X_data[j];
				}
				if (ABS(R_data[i][i]) > precision) {
					//检查R的对角元素是否零元。如果不是，可以直接求x。否则x不唯一，则令该行x元素为0 。
					X_data[i] = (QtB_data[i] - d) / R_data[i][i];
				}
				else {
					X_data[i] = 0;
				}
			}
		}

		//mat_print(Q);
		//mat_print(R);
		//mat_mul(Q, R, m);
		//mat_print(m);

		//vec_print(X);
		//vec_print(mX);

		mat_deallocate(&Q);
		mat_deallocate(&R);
		vec_deallocate(&QtB);

		return true;
	}
	else return false;
}



/// <summary>
/// 线性方程组最小二乘法解析求解。速度会比QR分解慢一些。
/// 如果m是横阵，这里求逆的时候，mtm的秩不够，会无法求逆。所以……横阵不要用解析法求最小二乘。
/// </summary>
/// <param name="m">矩阵m</param>
/// <param name="operation_buffer">运算的缓存矩阵。m.cols*3行，max(m.rows,m.cols)列。</param>
/// <param name="X">待求解的矢量X</param>
/// <param name="mX">已知的矢量mX</param>
/// <param name="precision">运算精度</param>
/// <returns>操作是否成功。若矩阵不匹配，或者矩阵不是满秩，运算失败，返回false</returns>
bool mat_LinearEquation_LinearLeastSquare_analytic(mat m, vec X, vec mX, real_t precision)
{
	if ((m.cols == X.size) && (m.rows == mX.size)) {  //先验证矩阵m和矢量X，mX的运算是否合法。

		//X有公式解： X = [invers(mt*m)]*mT*mX.      invers表示求逆。mt表示m的转置。
		//QR分解： 令m=q*r,mt=rt*qt，因为qt*q=I，所以X=[invers(mt*m)]*mT*mX=[invers(rt*r)]*rt*qt*mX。

		// (m.cols, m.cols);
		mat mtm;
		// (m.cols, m.cols);
		mat imtm;
		// (m.cols, m.rows);
		mat mt;
		// (m.cols, m.rows); //同mt，可以用mtm的内存。错峰使用
		mat imtmmt; 

		mat_allocate(&mtm, m.cols, m.cols);
		mat_allocate(&imtm, m.cols, m.cols);
		mat_allocate(&mt, m.cols, m.rows);
		mat_allocate(&imtmmt, m.cols, m.rows);


		mat_transpose_to(m, mt);
		mat_mul(mt, m, mtm);
 		// mat_print(mtm);

		mat_inversion(mtm, imtm, precision);  //如果m是横阵，这里求逆的时候，mtm的秩不够，会无法求逆。所以……横阵不要用解析法求最小二乘。

		mat_mul(imtm, mt, imtmmt);
		mat_mul_vec(imtmmt, mX, X);

		mat_deallocate(&mtm);
		mat_deallocate(&imtm);
		mat_deallocate(&mt);
		mat_deallocate(&imtmmt);

		return true;	
	}
	else return false;
}

/// <summary>
/// 矩阵乘法最小二乘法求解。m*X=mX，已知m和mX，求最佳的X。
/// 这里将矩阵拆分成列矢进行最小二乘求解。
/// </summary>
/// <param name="m">矩阵m</param>
/// <param name="operation_buffer">运算的缓存矩阵。m.rows*2+m.cols+3 行，max(m.rows+1,m.cols)列。</param>
/// <param name="X">待求解的矩阵X</param>
/// <param name="mX">已知的矩阵mX</param>
/// <param name="is_resetX">是否重置矢量X。如果重置,计算X前会初始化为零矢量</param>
/// <param name="precision">运算精度</param>
/// <returns>操作是否成功。若矩阵不匹配，或者矩阵不是满秩，运算失败，返回false</returns>
bool mat_LinearLeastSquare(mat m, mat X, mat mX, bool is_resetX, real_t precision)
{
	if ((m.cols == X.rows) && (m.rows == mX.rows) && (X.cols == mX.cols)) {  //先验证矩阵m和X，mX的运算是否合法。

		//operation_buffer: m.rows*2+m.cols+3 行，max(m.rows+1,m.cols)列。

		//m.cols
		vec vX;
		//m.rows
		vec vB;
		vec_allocate(&vX, m.cols);
		vec_allocate(&vB, m.rows);

		for (int i = 0; i < X.cols; i++) {
			//初始化vB
			for (int j = 0; j < vB.size; j++) {
				vB.data[j] = mX.data[j][i];
			}
			for (int j = 0; j < vX.size; j++) {
				vX.data[j] = X.data[j][i];
			}
			//printfs(
			mat_LinearEquation_LinearLeastSquare(m, vX, vB, is_resetX, precision);

			//mat_print(m);
			//vec_print(vB);
			//vec_print(vX);

			for (int j = 0; j < vX.size; j++) {
				X.data[j][i] = vX.data[j];
			}
		}
		vec_deallocate(&vX);
		vec_deallocate(&vB);

		return true;
	}
	else return false;
}


///////////  特征值分解        /////////////////

/// <summary>
/// 使用Householder矩阵将矩阵m化成Hessenberg形式。特征值分解的辅助函数。                                                                                                                                                                                                                                                                                                                                                       
/// </summary>
/// <param name="m">原矩阵=H*A*HT</param>
/// <param name="H">分解的正交矩阵</param>
/// <param name="A">分解的对三角矩阵，如果m不对称，则是上Hessenberg矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool EVD_Hessenberg_Householder(mat m,  mat H, mat A, real_t precision) //m矩阵=H*A*HT.
{
	//if (mat_isSymmetric(m,precision) == false) {
	//	return false;
	//}

	if ((m.rows == H.cols) && (m.rows == H.rows) &&
		(m.cols == A.cols) && (m.rows == A.rows))
	{
		//要cols个豪斯霍尔德向量，第一个向量rows个元素，第二个向量rows-1个元素.......
		//第i个向量rows-(i-1)个元素，最后一个向量rows-(cols-1-1)个元素。
		// 每个豪斯霍尔德向量需要一个beta参数。

		mat operation_buffer;
		vec beta_buffer;
		real_t** v;
		real_t beta[1];
		mat_allocate(&operation_buffer, m.cols, m.rows);
		vec_allocate(&beta_buffer, m.cols);

		v = operation_buffer.data;
		//内存分配完成。

		mat_copyto(m, A);//把该矩阵复制给A

		//下面开始变换
		int n = m.cols - 2;//选行列数比较小的值。

		for (int j = 0; j <= n; j++) {
			mat_Householder_left(A, j, j + 1, v[j], beta, false, precision);
			beta_buffer.data[j] = beta[0];
			mat_Householder_right(A, j + 1, j, v[j], beta, true, precision);
		}

		//下面计算H，H=H(0)*H(1)*……*H(n-2)*I。

		//初始化H为单位矩阵
		mat_eye(H);
		//开始计算H
		for (int j = n; j >= 0; j--) {
			beta[0] = beta_buffer.data[j];
			mat_Householder_left(H, j + 1, j + 1, v[j], beta, true, precision);
		}

		mat_deallocate(&operation_buffer);
		vec_deallocate(&beta_buffer);

		return true;
	}
	else return false;
}


/// <summary>
/// 使用Householder矩阵将矩阵m化成Hessenberg形式。特征值分解的辅助函数。                                                                                                                                                                                                                                                                                                                                                       
/// </summary>
/// <param name="m">原矩阵=H*A*HT</param>
/// <param name="A">分解的对三角矩阵，如果m不对称，则是上Hessenberg矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool EVD_Hessenberg_Householder_without_H(mat m, mat A, real_t precision) //m矩阵=H*A*HT.
{
	if ( (m.cols == A.cols) && (m.rows == A.rows))
	{
		vec operation_buffer;
		real_t* v;
		real_t beta[1];
		vec_allocate(&operation_buffer, m.rows);
		v = operation_buffer.data;
		//内存分配完成。

		mat_copyto(m, A);//把该矩阵复制给A

		//下面开始变换
		int n = m.cols - 2;//选行列数比较小的值。

		for (int j = 0; j <= n; j++) {

			mat_Householder_left(A, j, j + 1, v, beta, false, precision);
			mat_Householder_right(A, j + 1, j, v, beta, true, precision);
		}
		vec_deallocate(&operation_buffer);
		return true;
	}
	else return false;
}



/// <summary>
/// 快速排序(从小到大)。
/// </summary>
/// <param name="arr">需要排序的数组</param>
/// <param name="order">记录排序数组的序号，需要事先初始化</param>
/// <param name="begin">需要排序的区间左边界</param>
/// <param name="end">需要排序的区间的右边界，左右边界都取得到</param>
void quickSortUp(real_t* arr, int* order, int begin, int end)
{
	//如果区间不只一个数
	if (begin < end)
	{
		real_t temp = arr[begin];		//将区间的第一个数作为基准数,基准数先备份在temp
		int torder = order[begin];
		int i = begin;					//从左到右进行查找时的“指针”，指示当前左位置
		int j = end;						 //从右到左进行查找时的“指针”，指示当前右位置
		while (i < j)//这个循环内部的变量交换设计得很精巧。当这个大循环完成的时候，数据成功分成左右两份，右边的数据普遍比左边的大。
		{
			while (i < j && arr[j] > temp)  //当右边的数大于基准数时，略过，继续向左查找
				j--;
			arr[i] = arr[j];		//将右边小于等于基准元素的数填入左边相应位置,第一次i是begin的位置，它的数据已经备份在temp了。
			order[i] = order[j];
			while (i < j && arr[i] <= temp)
				i++;
			arr[j] = arr[i];      ////将左边大于基准元素的数填入右边相应位置。 j的数据刚刚传给这次循环前的i了	
			order[j] = order[i];
		}
		arr[i] = temp;		//将基准元素填入相应位置。此时的i即为基准元素的位置
		order[i] = torder;
		//对基准元素的左边子区间进行相似的快速排序
		quickSortUp(arr, order, begin, i - 1);
		//对基准元素的右边子区间进行相似的快速排序
		quickSortUp(arr, order, i + 1, end);
	}
	//如果区间只有一个数，则返回
	else
		return;
}

/// <summary>
/// 实现对称矩阵的本征值和本征向量的分解。 m = X*D*X^(T) 
/// </summary>
/// <param name="m">待分解的对称矩阵</param>
/// <param name="operation_buffer">运算的缓存矩阵，m.rows*3+1行，m.rows+1列。</param>
/// <param name="X">X矩阵的列向量是本征向量</param>
/// <param name="D">对角矩阵，对角元素是本征值。</param>
/// <param name="precision">运算精度</param>
/// <returns>迭代计算的次数。如果不是对称矩阵，计算失败，返回-1</returns>
int mat_EVD_Symmetric(mat m,  mat X, mat D, real_t precision)
{
	if (mat_isSymmetric(m, precision) == false) {
		return -1;
	}

	if (EVD_Hessenberg_Householder(m, X, D, precision)) {		
		//EVD_Hessenberg_Householder函数里有检查矩阵的尺寸了。
		
		//现在X是单位正交矩阵，D是三对角矩阵。
		

		//带原点位移QR迭代。迭代2*n次。
		//QR分解需要的变量
		//Matrix H(this.rows, this.rows);
		//Matrix Q(this.rows, this.rows);
		//Matrix R(this.rows, this.rows);

		mat H, Q, R;  //下面要给它们分配内存。
		ivec order_buffer;
		mat_allocate(&H, m.rows, m.rows);
		mat_allocate(&Q, m.rows, m.rows);
		mat_allocate(&R, m.rows, m.rows);
		ivec_allocate(&order_buffer, m.rows);

		real_t** H_data = H.data;
		real_t** Q_data = Q.data;
		real_t** R_data = R.data;
		real_t** X_data = X.data;
		real_t** D_data = D.data;

		real_t u = 0;//平移量
		real_t d = 0;//用于计算更好的平移量
		int i = 0, j = 0;
		mat_copyto(X, H);//X,H,相互更新。

		int times = 0;

		int sm = 0, sn = 0; //记录qr分解的子矩阵的位置。从sm行列到sn行列。
		bool is_diag = false;
		bool flag = false;
		bool isget_mn = false;
		int sub_times;  //记录子循环的次数。
		while (is_diag == false) //下面分块进行QR分解迭代。
		{
			//设计思路:设初始flag为0，依次检测斜对角元素，为零flag置1，否则置0。检测上升沿和下降沿。
			//上升沿设置sm，下降沿设置sn。sm，sn是QR迭代的子三对角矩阵区域。设置了mn即可开始迭代。
			//如果无下降沿，则sn=m.cols-1。如无上升沿，则斜对角线都是零元，可结束迭代。

			//初始化标志位
			flag = false;
			isget_mn = false;
			sub_times = 0;
			for (i = 1; i < m.cols; i++) {
				if (ABS(D_data[i - 1][i]) > precision) {
					if (flag == false) {	//检测上升沿。
						sm = i - 1;
					}
					if (i == m.cols - 1) {
						sn = i;
						isget_mn = true;
					}
					flag = true;
				}
				else {
					if (flag == true) {		//检测下降沿
						sn = i - 1;
						isget_mn = true;
					}
					flag = false;
				}

				if (isget_mn) {
					//QR分解迭代，如果非对角元素绝对值大于precision，进行平移迭代。
					//计算平移量。该计算公式见《数值线性代数，徐树方》P208	       //另一个简单的平移量：u = D_data[i][i]; 效果之比这个差一点
					d = (D_data[sn - 1][sn - 1] - D_data[sn][sn]) / 2;
					if (d > 0) {
						u = D_data[sn][sn] + d - sqrtf(d * d + D_data[sn - 1][sn] * D_data[sn - 1][sn]);
					}
					else {
						u = D_data[sn][sn] + d + sqrtf(d * d + D_data[sn - 1][sn] * D_data[sn - 1][sn]);
					}
					//平移
					for (j = sm; j <= sn; j++) {
						D_data[j][j] = D_data[j][j] - u;
					}
					//Q.eye(sm,sn);	
					{
						for (int ii = sm; ii <= sn; ii++) {
							for (int ij = sm; ij < i; ij++) {
								Q_data[ij][ii] = 0;
								Q_data[ii][ij] = 0;
							}
							Q_data[ii][ii] = 1;
						}
					}				
					//  D在EVD_subfunction3_QR_Givens里被改变
					// D.EVD_subfunction3_QR_Givens( Q, sm, sn, precision);	//D是Hessenberg矩阵，下三角区域几乎全是零。这里用Givens变换更快
					{
						//EVD_subfunction3_QR_Givens 函数实现
						real_t** this_data = D_data;
						real_t** G_data = Q_data;
						int n0 = sm;
						int n1 = sn;
						//
					//	int i, j;
						int i2, j2;
						real_t sin, cos; //保存吉文斯矩阵参数
						real_t d0, d1, d2;	 //d1,d2用来保存矩阵在进行吉文斯变换时的中间参数

						//直接改变this的矩阵。this   <<---  G * this
						//开始Givens变换
						{
							for (j = n0; j <= n1; j++) {		//逐列消元
								for (i = j + 1; i <= n1; i++) {  //消下三角区域元素
									if (ABS(this_data[i][j]) >= precision) {  //只有下三角元素不为零才需要吉文斯变换
										//计算吉文斯矩阵参数
										d0 = sqrtf(this_data[j][j] * this_data[j][j] + this_data[i][j] * this_data[i][j]);
										cos = this_data[j][j] / d0;    //G= |cos,  sin|  ,G是单位正交对称矩阵，G^(-1)=G
										sin = this_data[i][j] / d0;	//		|sin,-cos|
										//更新A：A左乘G(i).  A=G(i)*A
										for (j2 = j; (j2 <= j + 2) && (j2 <= n1); j2++) {    //从第j列开始计算。因为0--(j-1)列下三角都是0了，不用去计算。
											d1 = this_data[j][j2];			//因为这是对三角矩阵，A右上角有大量零元区	
											d2 = this_data[i][j2];
											this_data[j][j2] = d1 * cos + d2 * sin;
											this_data[i][j2] = d1 * sin - d2 * cos;
										}
										//更新G：G右乘G(i).  G=G*G(i)  
										for (i2 = n0; i2 <= i; i2++) {		//有优化空间。可以在草稿纸上记录零元区域的变换情况，从而避免不必要的零元计算。
											d1 = G_data[i2][j];						//零元区域变化是有规律的。这里已经避免了零元区域的计算了。
											d2 = G_data[i2][i];
											G_data[i2][j] = d1 * cos + d2 * sin;
											G_data[i2][i] = d1 * sin - d2 * cos;
										}
									}
								}
							}
						}
					}


					//	X.mul(Q, H);    // X和H是相互更新的。
					for ( i = 0; i < X.rows; i++) {	//R是sm行列到sn行列的上三角，其他地方是单位或者0矩阵。
						for ( j = sm; j <= sn; j++) {
							d = 0;
							for (int k = sm; k <= sn; k++) {
								d += (X_data[i][k] * Q_data[k][j]);
							}
							H_data[i][j] = d;	//访问B_data[i][j]要先访问B,再移到data，再到data[i]，最后到data[i][j]，寻址效率比较低。所以用中间变量来提高速度。
						}
					}

					//  D.mul(Q, R);   //有大量零元，可以专门设计乘法以减少计算量。这里D与R相互更新。						
					{//	D.copyto(R);   //注意，在吉文斯变换中，Q只是局部地变换。	
						for ( i = 1; i < D.rows - 2; i++) {
							R_data[i][i - 1] = D_data[i][i - 1];
							R_data[i][i] = D_data[i][i];
							R_data[i][i + 1] = D_data[i][i + 1];
							R_data[i][i + 2] = D_data[i][i + 2];
						}
						//这时i = D.rows - 2，R_data[i][i + 2] = D_data[i][i + 2]; 的语句会内存溢出
						R_data[i][i - 1] = D_data[i][i - 1];  
						R_data[i][i] = D_data[i][i];
						R_data[i][i + 1] = D_data[i][i + 1];

						R_data[0][0] = D_data[0][0];
						R_data[0][1] = D_data[0][1];
						R_data[0][2] = D_data[0][2];
						R_data[R.rows - 1][R.rows - 2] = D_data[R.rows - 1][R.rows - 2];
						R_data[R.rows - 1][R.rows - 1] = D_data[R.rows - 1][R.rows - 1];
					}
					for ( i = 0; i < R.rows; i++) {	//实现乘法R.mul(Q, D);。
						for ( j = sm; j <= sn; j++) {
							d = 0;
							for (int k = (i > sm) ? (i) : (sm); k <= sn; k++) {
								d += (D_data[i][k] * Q_data[k][j]);
							}
							R_data[i][j] = d;	//访问B_data[i][j]要先访问B,再移到data，再到data[i]，最后到data[i][j]，寻址效率比较低。所以用中间变量来提高速度。
						}
					}

					//平移回去。注意，D与R相互更新。现在是要平移R。
					for (j = sm; j <= sn; j++) {
						R_data[j][j] = R_data[j][j] + u;
					}

					//这里一次循环相当于做了两次QR分解迭代

					//计算平移量
					d = (R_data[sn - 1][sn - 1] - R_data[sn][sn]) / 2;
					if (d > 0) {
						u = R_data[sn][sn] + d - sqrtf(d * d + R_data[sn - 1][sn] * R_data[sn - 1][sn]);
					}
					else {
						u = R_data[sn][sn] + d + sqrtf(d * d + R_data[sn - 1][sn] * R_data[sn - 1][sn]);
					}
					//平移
					for (j = sm; j <= sn; j++) {
						R_data[j][j] = R_data[j][j] - u;
					}
					//R.EVD_subfunction3_QR_Givens(Q, sm, sn, precision);

					//Q.eye(sm,sn);	
					{
						for ( i = sm; i <= sn; i++) {
							for ( j = sm; j < i; j++) {
								Q_data[j][i] = 0;
								Q_data[i][j] = 0;
							}
							Q_data[i][i] = 1;
						}
					}
					// R.EVD_subfunction3_QR_Givens( Q, sm, sn, precision);	//D是Hessenberg矩阵，下三角区域几乎全是零。这里用Givens变换更快
					{
						//EVD_subfunction3_QR_Givens 函数实现
						real_t** this_data = R_data;
						real_t** G_data = Q_data;
						int n0 = sm;
						int n1 = sn;
						//
						//int i, j;
						int i2, j2;
						real_t sin, cos; //保存吉文斯矩阵参数
						real_t d0, d1, d2;	 //d1,d2用来保存矩阵在进行吉文斯变换时的中间参数

						//直接改变this的矩阵。this   <<---  G * this
						//开始Givens变换
						{
							for (j = n0; j <= n1; j++) {		//逐列消元
								for (i = j + 1; i <= n1; i++) {  //消下三角区域元素
									if (ABS(this_data[i][j]) >= precision) {  //只有下三角元素不为零才需要吉文斯变换
										//计算吉文斯矩阵参数
										d0 = sqrtf(this_data[j][j] * this_data[j][j] + this_data[i][j] * this_data[i][j]);
										cos = this_data[j][j] / d0;    //G= |cos,  sin|  ,G是单位正交对称矩阵，G^(-1)=G
										sin = this_data[i][j] / d0;	//		|sin,-cos|
										//更新A：A左乘G(i).  A=G(i)*A
										for (j2 = j; (j2 <= j + 2) && (j2 <= n1); j2++) {    //从第j列开始计算。因为0--(j-1)列下三角都是0了，不用去计算。
											d1 = this_data[j][j2];			//因为这是对三角矩阵，A右上角有大量零元区	
											d2 = this_data[i][j2];
											this_data[j][j2] = d1 * cos + d2 * sin;
											this_data[i][j2] = d1 * sin - d2 * cos;
										}
										//更新G：G右乘G(i).  G=G*G(i)  
										for (i2 = n0; i2 <= i; i2++) {		//有优化空间。可以在草稿纸上记录零元区域的变换情况，从而避免不必要的零元计算。
											d1 = G_data[i2][j];						//零元区域变化是有规律的。这里已经避免了零元区域的计算了。
											d2 = G_data[i2][i];
											G_data[i2][j] = d1 * cos + d2 * sin;
											G_data[i2][i] = d1 * sin - d2 * cos;
										}
									}
								}
							}
						}
					}


					//H.mul(Q, X);
					for ( i = 0; i < H.rows; i++) {	//R是sm行列到sn行列的上三角，其他地方是单位或者0矩阵。
						for ( j = sm; j <= sn; j++) {
							d = 0;
							for (int k = sm; k <= sn; k++) {
								d += (H_data[i][k] * Q_data[k][j]);
							}
							X_data[i][j] = d;	//访问B_data[i][j]要先访问B,再移到data，再到data[i]，最后到data[i][j]，寻址效率比较低。所以用中间变量来提高速度。
						}
					}
					//	R.mul(Q, D); //上三角矩阵乘方阵，可以用专门的乘法减少计算量.				
					{//	R.copyto(D);
						for ( i = 1; i < R.rows - 2; i++) {
							D_data[i][i - 1] = R_data[i][i - 1];
							D_data[i][i] = R_data[i][i];
							D_data[i][i + 1] = R_data[i][i + 1];
							D_data[i][i + 2] = R_data[i][i + 2];
						}
						//这时i = D.rows - 2，D_data[i][i + 2] = R_data[i][i + 2]; 的语句会内存溢出
						D_data[i][i - 1] = R_data[i][i - 1];
						D_data[i][i] = R_data[i][i];
						D_data[i][i + 1] = R_data[i][i + 1];

						D_data[0][0] = R_data[0][0];
						D_data[0][1] = R_data[0][1];
						D_data[0][2] = R_data[0][2];
						D_data[R.rows - 1][R.rows - 2] = R_data[R.rows - 1][R.rows - 2];
						D_data[R.rows - 1][R.rows - 1] = R_data[R.rows - 1][R.rows - 1];
					}
					for ( i = 0; i < R.rows; i++) {
						for ( j = sm; j <= sn; j++) {
							d = 0;
							for (int k = (i > sm) ? (i) : (sm); k <= sn; k++) {
								d += (R_data[i][k] * Q_data[k][j]);
							}
							D_data[i][j] = d;	//访问B_data[i][j]要先访问B,再移到data，再到data[i]，最后到data[i][j]，寻址效率比较低。所以用中间变量来提高速度。
						}
					}
					//平移回去
					for (j = sm; j <= sn; j++) {
						D_data[j][j] = D_data[j][j] + u;
					}

					//迭代一次后就重置标记
					isget_mn = false;
					sub_times = sub_times + 2;
					//	printfs(" sub_times =  " << sub_times << "  。 \r\n";
				}
			}
			//如果扫描了所有斜对角线元素，都没有确定一个mn，则说明全是零元，无上升沿。分解成对角矩阵了。
			if (sub_times == 0) {
				is_diag = true;
			}
			times = times + 2;
			//	printfs("    这是第     " << times << "   次迭代   \r\n";

		}

		//给特征值排序。使其从大到小排列。
		{
			int length=m.rows; //排序数组的长度。

			//次序。排序之前值是0--length-1。order跟待排序的数组一同变换。
			//以记录从小到大排序的第n个数字在排序前的序号。

			int* order=order_buffer.data; 
			real_t* lamda = R.data[0];	    //矩阵R后面用不到了。就拿内存来辅助排序吧。其实还是operation_buffer的内存
			real_t* abslamda = R.data[1]; 
			//这里特征值分解的最小矩阵 是2*2的矩阵。 不然这里内存会溢出。

			for ( i = 0; i < length; i++) {
				lamda[i] = D_data[i][i];
				abslamda[i] = ABS(D_data[i][i]);
				order[i] = i;
			}
			//下面对abslamda进行排序.
			quickSortUp(abslamda, order, 0, length - 1);

			//现在根据特征值大小从大到小组织特征向量矩阵。
			for ( j = 0; j < m.cols; j++) {
				for ( i = 0; i < m.rows; i++) {
					H_data[i][j] = X_data[i][order[m.cols - 1 - j]];
				}
			}
			mat_copyfrom(X, H);
			mat_eye(D);
			for ( i = 0; i < m.cols; i++) {
				D_data[i][i] = lamda[order[m.cols - 1 - i]];
			}
		}

		mat_deallocate(&H);
		mat_deallocate(&Q);
		mat_deallocate(&R);
		ivec_deallocate(&order_buffer);

		return times;
	}
	else return (-1);
}



///////////  奇异值分解        /////////////////

/// <summary>
/// 实现矩阵的二对角化。奇异值分解的子函数。
/// </summary>
/// <param name="m">原矩阵m= U * B * VT</param>
/// <param name="operation_buffer">运算缓存矩阵。行数m.col*2；列数max（m.col，m.row）。奇异值分解已经默认行数大于列数了</param>
/// <param name="U">单位正交矩阵。m.rows行列的方阵</param>
/// <param name="B">是对二角矩阵。大小与m一致</param>
/// <param name="VT">单位正交矩阵。m.cols行列的方阵</param>
/// <param name="precision">运算精度</param>
/// <returns>操作是否成功。若矩阵不匹配，返回false</returns>
bool SVD_subfunction1(mat m, mat U, mat B, mat VT, real_t precision)
{
	//this矩阵m*n(m行n列),U是m*m，B是m*n,V是n*n。
	if ((m.rows == U.rows) && (m.rows == U.cols) &&
		(m.rows == B.rows) && (m.cols == B.cols) &&
		(m.cols == VT.rows) && (m.cols == VT.cols)) {
		//思路：用《数值线性代数，徐树方》P233	的方法。 默认行数大于列数。如果列数大于行数，则退出程序。
		if (m.rows < m.cols) {
			return false;
		}
		// operation_buffer 列数：max（m.col，m.row），行数：m.col*2行
		int n = B.cols;
		real_t beta[1];
		mat u_buffer;
		mat v_buffer;
		vec betau_buffer;
		vec betav_buffer;
		mat_allocate(&u_buffer, m.cols, m.rows);
		mat_allocate(&v_buffer, m.cols-2, m.cols);
		vec_allocate(&betau_buffer, n);
		vec_allocate(&betav_buffer, n);

		//先分配内存。
		//矩阵U要n个豪斯霍尔德向量，第一个向量m.rows个元素，第二个向量m.rows-1个元素.......
		//第i个向量m.rows-i+1个元素，最后一个向量m.rows-n+1个元素。
		real_t** u = u_buffer.data ;       //col行，row列
		real_t* u_beta = betau_buffer.data; //1行，col列
		//矩阵VT要n-2个豪斯霍尔德向量，第一个向量m.cols-1个元素，第二个向量m.cols-2个元素.......
		//第i个向量m.cols-i个元素，最后一个向量m.cols-n+2个元素。
		real_t** v = v_buffer.data; // col-2行，col列
		real_t* v_beta = betav_buffer.data; // 1行，col列。

		int i, j;
		//把该矩阵复制给B
		//mat_copyto(m, B);  //在外面已经复制过了

		//下面开始变换
		{
			for (j = 0; j < n - 2; j++) {		//行数大于列数，正常情况
				i = j + 1;

				mat_Householder_left(B, j, j, u[j], beta, false, precision);
				u_beta[j] = beta[0];
				mat_Householder_right(B, i, j, v[j], beta, false, precision);
				v_beta[j] = beta[0];
			}

			j = n - 2;
			mat_Householder_left(B, j, j, u[j], beta, false, precision);
			u_beta[j] = beta[0];

			j = n - 1;
			mat_Householder_left(B, j, j, u[j], beta, false, precision);
			u_beta[j] = beta[0];
		}

		//发现大矩阵的奇异值分解中，矩阵二对角化需要消耗相当多的时间，这里可以输出字符串来显示一下进度。

			//下面计算U和VT，H=H(0)*H(1)*……*H(n-2)*I。
		mat_eye(U);//初始化
		mat_eye(VT);//初始化
		{
			for (j = n - 1; j >= 0; j--) {
				beta[0] = u_beta[j];
				mat_Householder_left(U, j, j, u[j], beta, true, precision);
			}
			for (j = n - 2; j > 0; j--) {
				beta[0] = v_beta[j - 1];
				mat_Householder_right(VT, j, j, v[j - 1], beta, true, precision);
			}
		}

		mat_deallocate(&u_buffer);
		mat_deallocate(&v_buffer);
		vec_deallocate(&betau_buffer);
		vec_deallocate(&betav_buffer);

		return true;
	}
	else return false;
}

/// <summary>
/// 实现《数值线性代数，徐树方》P234的QR迭代方法。
/// 调用一次该函数，对B的m到n的子矩阵进行一次迭代。 B <<-- P*B*Q
/// </summary>
/// <param name="P">是SVD_subfunction1的U</param>
/// <param name="B">是SVD_subfunction1的上二对角矩阵</param>
/// <param name="Q">是SVD_subfunction1的VT</param>
/// <param name="m">子矩阵起始行列</param>
/// <param name="n">子矩阵结束行列</param>
/// <param name="precision">运算精度</param>
void SVD_subfunction2(mat P, mat B, mat Q, int m, int n, real_t precision)
{
	real_t** P_data = P.data;
	real_t** B_data = B.data;
	real_t** Q_data = Q.data;

	int i, j;
	real_t alpha, beta, deta, u, d;
	real_t x, y;
	real_t sin, cos;//吉文斯变换参数。
	//先确定平移量u，（这里的平移量其实并不需要真的平移，只是辅助计算最佳的第一个sin和cos）
	alpha = B_data[n][n] * B_data[n][n] + B_data[n - 1][n] * B_data[n - 1][n];
	if (n - 2 >= 0) {
		deta = (B_data[n - 1][n - 1] * B_data[n - 1][n - 1] + B_data[n - 2][n - 1] * B_data[n - 2][n - 1] - alpha) / 2;
	}
	else {
		deta = (B_data[n - 1][n - 1] * B_data[n - 1][n - 1] - alpha) / 2;
	}
	beta = B_data[n - 1][n] * B_data[n - 1][n - 1];
	if (deta > 0) {
		u = alpha - beta * beta / (deta + sqrtf(deta * deta + beta * beta));
	}
	else {
		u = alpha - beta * beta / (deta - sqrtf(deta * deta + beta * beta));
	}
	x = B_data[m][m] * B_data[m][m] - u;
	y = B_data[m][m] * B_data[m][m + 1];
	d = sqrtf(x * x + y * y);
	sin = y / d;		//	G = |cos   sin|
	cos = x / d;		//      |sin  -cos|

	//第一轮QR迭代。	//  B <<-- B*G
	x = B_data[m][m];
	y = B_data[m][m + 1];
	B_data[m][m] = x * cos + y * sin;
	B_data[m][m + 1] = x * sin - y * cos;
	x = B_data[m + 1][m];
	y = B_data[m + 1][m + 1];
	B_data[m + 1][m] = x * cos + y * sin;
	B_data[m + 1][m + 1] = x * sin - y * cos;
	//第一次更新Q
	for (j = 0; j < Q.cols; j++) {
		x = Q_data[m][j];
		y = Q_data[m + 1][j];
		Q_data[m][j] = x * cos + y * sin;
		Q_data[m + 1][j] = x * sin - y * cos;
	}

	//进入QR迭代。去除第一轮和最后一轮，中间的QR迭代都很有规律。
	for (i = m; i <= n - 2; i++) {
		//  B <<--- G*B
		x = B_data[i][i];
		y = B_data[i + 1][i];
		d = sqrtf(x * x + y * y);
		sin = y / d;		//	G = |cos   sin|
		cos = x / d;		//      |sin  -cos|
		for (j = i; j < i + 3; j++) {
			x = B_data[i][j];
			y = B_data[i + 1][j];
			B_data[i][j] = x * cos + y * sin;
			B_data[i + 1][j] = x * sin - y * cos;
		}
		//更新P。P <<-- I*G2*G4*G6……  
		for (j = 0; j < P.rows; j++) {
			x = P_data[j][i];
			y = P_data[j][i + 1];
			P_data[j][i] = x * cos + y * sin;
			P_data[j][i + 1] = x * sin - y * cos;
		}


		//  B <<--- B*G
		x = B_data[i][i + 1];
		y = B_data[i][i + 2];
		d = sqrtf(x * x + y * y);
		sin = y / d;		//	G = |cos   sin|
		cos = x / d;		//      |sin  -cos|
		for (j = i; j < i + 3; j++) {
			x = B_data[j][i + 1];
			y = B_data[j][i + 2];
			B_data[j][i + 1] = x * cos + y * sin;
			B_data[j][i + 2] = x * sin - y * cos;
		}
		//更新Q。Q <<--  ……G5*G3*G1*I  
		for (j = 0; j < Q.cols; j++) {
			x = Q_data[i + 1][j];
			y = Q_data[i + 2][j];
			Q_data[i + 1][j] = x * cos + y * sin;
			Q_data[i + 2][j] = x * sin - y * cos;
		}
	}

	//最后一轮QR迭代。	//  B <<--- G*B
	x = B_data[n - 1][n - 1];
	y = B_data[n][n - 1];
	d = sqrtf(x * x + y * y);
	sin = y / d;		//	G = |cos   sin|
	cos = x / d;		//      |sin  -cos|
	x = B_data[n - 1][n - 1];
	y = B_data[n][n - 1];
	B_data[n - 1][n - 1] = x * cos + y * sin;
	B_data[n][n - 1] = x * sin - y * cos;
	x = B_data[n - 1][n];
	y = B_data[n][n];
	B_data[n - 1][n] = x * cos + y * sin;
	B_data[n][n] = x * sin - y * cos;
	//更新P。P <<-- I*G2*G4*G6……  
	for (j = 0; j < P.rows; j++) {
		x = P_data[j][n - 1];
		y = P_data[j][n];
		P_data[j][n - 1] = x * cos + y * sin;
		P_data[j][n] = x * sin - y * cos;
	}

}

/// <summary>
/// 实现《数值线性代数，徐树方》P237的QR变换方法。
/// 调用一次该函数，将B的左上方块的对角零元转换成零元行。 B <<--  U*B。  默认行数大于列数
/// </summary>
/// <param name="U">来记录吉文斯变换的矩阵</param>
/// <param name="B">上二对角矩阵</param>
/// <param name="precision">运算精度</param>
void SVD_subfunction3( mat U, mat B, real_t precision)
{
	real_t** U_data = U.data;
	real_t** B_data = B.data;

	int i, j;
	int m=-1;   //记录零元行。
	real_t x, y, d;
	real_t sin, cos;	//吉文斯变换参数。
	bool isget_m = false;
	//先检测对角元，看有无零元。并检测该零元是否在零元行。
	for (i = 0; i < B.cols - 1; i++) {
		if ((ABS(B_data[i][i]) < precision) && (ABS(B_data[i][i + 1]) > precision)) {
			m = i;
			isget_m = true;
			break;
		}
	}
	if (isget_m) {
		for (i = m + 1; (i < B.cols - 1) && (ABS(B_data[i][i + 1]) > precision); i++) {
			x = B_data[m][i];
			y = B_data[i][i];
			d = sqrtf(x * x + y * y);
			cos = -y / d;		//	G = |cos   sin|		.这里情况有变，是要使x为0.所以cos=-sin，sin=cos。
			sin = x / d;		//      |sin  -cos|
			{
				x = B_data[m][i];
				y = B_data[i][i];
				B_data[m][i] = x * cos + y * sin;
				B_data[i][i] = x * sin - y * cos;
				x = B_data[m][i + 1];
				y = B_data[i][i + 1];
				B_data[m][i + 1] = x * cos + y * sin;
				B_data[i][i + 1] = x * sin - y * cos;
			}
			//更新U。U <<-- I*G2*G4*G6……  
			for (j = 0; j < U.rows; j++) {
				x = U_data[j][m];
				y = U_data[j][i];
				U_data[j][m] = x * cos + y * sin;
				U_data[j][i] = x * sin - y * cos;
			}
		}

		//最后一列变换。i已经在最后一次循环里加1了。
		x = B_data[m][i];
		y = B_data[i][i];
		d = sqrtf(x * x + y * y);
		cos = -y / d;		//	G = |cos   sin|
		sin = x / d;		//      |sin  -cos|
		{
			x = B_data[m][i];
			y = B_data[i][i];
			B_data[m][i] = x * cos + y * sin;
			B_data[i][i] = x * sin - y * cos;
		}
		//更新U。U <<-- I*G2*G4*G6……  
		for (j = 0; j < U.rows; j++) {
			x = U_data[j][m];
			y = U_data[j][i];
			U_data[j][m] = x * cos + y * sin;
			U_data[j][i] = x * sin - y * cos;
		}
	}

}

/// <summary>
/// 实现奇异值分解
/// </summary>
/// <param name="m">m矩阵 = m_U * m_B * m_VT .</param>
/// <param name="operation_buffer">运算缓存矩阵。令：max=max(row,col)，min=min(row,col)，缓存矩阵要max*2+min*3行，max列</param>
/// <param name="m_U">m.rows*m.rows的单位正交矩阵，列矢量之间单位正交</param>
/// <param name="m_B">对角矩阵。大小和矩阵m一样，对角元素是奇异值。</param>
/// <param name="m_VT">m.cols*m.cols的单位正交矩阵，行矢量之间单位正交</param>
/// <param name="precision">运算精度</param>
/// <returns>迭代次数</returns>
int mat_SVD(mat m_m,  mat m_U, mat m_B, mat m_VT, real_t precision)
{
	//行数大于列数。否则就对this矩阵的转置做奇异值分解。最后奇异值排序时转置回来。
	int times = 0;
	int row, col;
	if (m_m.rows < m_m.cols) {
		row = m_m.cols;
		col = m_m.rows;
	}
	else {
		row = m_m.rows;
		col = m_m.cols;
	}

	// 计算 operation_buffer内存大小：令： max=max(row,col)，min=min(row,col)。
	// 对二角化要min*2行，max列。复制矩阵需要max*2+min行。共max*2+min*3行

	//要准备this矩阵，U矩阵，VT矩阵的备份，排列奇异值的时候会用到。 备份是行数大于列数的。
	mat U; U.rows = row; U.cols = row; 
	mat B; B.rows = row; B.cols = col; 
	mat VT; VT.rows = col; VT.cols = col;
	ivec order_buffer;
	vec lamda_buffer;
	vec abslamda_buffer;

	mat_allocate(&U, row, row);
	mat_allocate(&B, row, col);
	mat_allocate(&VT, col, col);
	ivec_allocate(&order_buffer, m_m.cols);  //行数大于列数，奇异值没那么多
	vec_allocate(&lamda_buffer, m_m.cols);
	vec_allocate(&abslamda_buffer, m_m.cols);

	real_t** B_data = B.data;
	if (m_m.rows < m_m.cols) {
		mat_transpose_to(m_m, B);
	}
	else {
		mat_copyto(m_m, B);
	}

	//this矩阵m*n(m行n列)，U是m*m，S是m*n,V是n*n。 行数大于列数。
	if (SVD_subfunction1(B, U, B, VT, precision)) {
		//printfs("    成功实现二对角化 ,打印对二角矩阵      \r\n";
		//mat_print(B);
		//mat_print(U);
		//mat_print(vT);

		//	printfs("    成功实现二对角化       \r\n";
		int i = 0, j = 0;
		int m = 0, n = 0; //记录qr分解的子矩阵的位置。从m行列到n行列。
		bool is_diag = false;
		bool flag = false;
		bool isget_mn = false;
		int sub_times;
		while (is_diag == false) //下面分块进行QR分解迭代。
		{
			SVD_subfunction3(U, B, precision);  
			


			//先检查对角元有无零元。

			//设计思路:设初始flag为0，依次检测斜对角元素，为零flag置1，否则置0。检测上升沿和下降沿。
			//上升沿设置m，下降沿设置n。m，n是QR迭代的子三对角矩阵区域。设置了mn即可开始迭代。
			//如果无下降沿，则n=m.cols-1。如无上升沿，则斜对角线都是零元，可结束迭代。

			//初始化标志位
			flag = false; //是否检测到上升沿的标志位。
			isget_mn = false;
			sub_times = 0;
			//将二对角矩阵分解成多个分块的单位矩阵，可以减少计算量。
			for (i = 1; i < B.cols; i++) {
				if (ABS(B_data[i - 1][i]) > precision) {
					if (flag == false) {
						m = i - 1;
					}
					if (i == (B.cols - 1)) {  //如果i是次对角线最后一个元素，则n的值也可以知道。
						n = i;
						isget_mn = true;
					}
					flag = true;
				}
				else {
					if (flag == true) {		//检测下降沿
						n = i - 1;
						isget_mn = true;
						flag = false;
					}
				}
				if (isget_mn) {
					SVD_subfunction2(U, B, VT, m, n, precision);
					//迭代一次后就重置标记
					isget_mn = false;
					sub_times = sub_times + 2;
					//	printfs(" sub_times =  " << sub_times << "  。 \r\n";
				}
			}

			//如果扫描了所有斜对角线元素，都没有确定一个mn，则说明全是零元，无上升沿。分解成对角矩阵了。
			if (sub_times == 0) {
				is_diag = true;

				//printfs(" \r\n  print U B VT \r\n");
				//mat_print(U);
				//mat_print(B);
				//mat_print(VT);

			}
			times = times + 1;
			//	printfs("    这是第     " << times << "   次迭代   \r\n";
		}

		//下面排列奇异值。共col个奇异值

		//将对二角化的运算缓存拿来辅助排序。
		real_t* lamda = lamda_buffer.data; 
		real_t* abslamda = abslamda_buffer.data;
		int* order = order_buffer.data;
		for ( i = 0; i < col; i++) {
			lamda[i] = B_data[i][i];
			abslamda[i] = ABS(B_data[i][i]);
			order[i] = i;
		}
		quickSortUp(abslamda, order, 0, col - 1);
		//现在根据奇异值大小从大到小组织矩阵。
		mat_init(m_B, 0);
		for ( i = 0; i < col; i++) {
			m_B.data[i][i] = lamda[order[col - 1 - i]];
		}
		real_t** m_U_data = m_U.data;
		real_t** m_VT_data = m_VT.data;
		real_t** U_data = U.data;
		real_t** VT_data = VT.data;

		if (m_B.rows < m_B.cols) {	//如果列数大于行数，说明矩阵是转置后再分解的。U是m_VT的转置。VT是m_U的转置。记得转置
			for ( j = 0; j < col; j++) {
				for ( i = 0; i < VT.rows; i++) {
					m_U_data[i][j] = VT_data[order[col - 1 - j]][i];
				}
				for ( i = 0; i < U.rows; i++) {
					m_VT_data[j][i] = U_data[i][order[col - 1 - j]];
				}
			}
			for ( j = col; j < U.cols; j++) {
				for ( i = 0; i < U.rows; i++) {
					m_VT_data[j][i] = U_data[i][j];
				}
			}
		}
		else {
			for ( j = 0; j < col; j++) {
				for ( i = 0; i < U.rows; i++) {
					m_U_data[i][j] = U_data[i][order[col - 1 - j]];
				}
				for ( i = 0; i < VT.cols; i++) {
					m_VT_data[j][i] = VT_data[order[col - 1 - j]][i];
				}
			}
			for ( j = col; j < U.cols; j++) {
				for ( i = 0; i < U.rows; i++) {
					m_U_data[i][j] = U_data[i][j];
				}
			}
		}
	}
	else {
		times = -1;
	}

	mat_deallocate(&U);
	mat_deallocate(&B);
	mat_deallocate(&VT);
	ivec_deallocate(&order_buffer);
	vec_deallocate(&lamda_buffer);
	vec_deallocate(&abslamda_buffer);
	return times;

}


// 非对称矩阵的特征值分解 会出现复数的特征值。这个函数很有用，可以用来计算任意阶数方程的根。



/// <summary>
/// 豪斯霍尔德变换。
/// 变换后，A = H * A 。H是豪斯霍尔德矩阵，以豪斯霍尔德矢量ｖ和参数beta的形式出现。 H=I-(1/beta)v*v^(T)
/// 变换后的A的x1列，y1行以下到y2的元素是0元素，把x1列y1行开始的列向量变成(a,0,0,...)^(T)形式的向量。得到的豪斯霍尔德矩阵以向量形式保存在v中。
/// 注意事项：注意数组v的大小，防止越界。
/// 这里是前闭后开。即，前面取得到，后面取不到
/// </summary>
/// <param name="A">矩阵</param>
/// <param name="x1">矩阵豪斯霍尔德变换的起始列</param>
/// <param name="x2">矩阵豪斯霍尔德变换的终止列</param>
/// <param name="y1">矩阵豪斯霍尔德变换的起始行</param>
/// <param name="y2">矩阵豪斯霍尔德变换的终止行</param>
/// <param name="v">豪斯霍尔德矢量，大小是变化的行数量</param>
/// <param name="beta">豪斯霍尔德矢量的参数</param>
/// <param name="isH_ensured">是否已经确认豪斯霍尔德矢量，如果未确认，函数会先计算更新豪斯霍尔德矢量和参数</param>
/// <param name="precision">计算精度</param>
void mat_evd_subf_Householder_left(mat A, int x1,int x2, int y1, int y2, real_t* v, real_t* beta, bool isH_ensured, real_t precision)
{
	int i, j;
	real_t d;
	real_t** A_data = A.data;

	if (x2 > A.cols)x2 = A.cols;
	if (y2 > A.rows)y2 = A.rows;
	if (isH_ensured == false) {  //如果还没确定Householder向量，先计算Householder向量。
		d = 0;
		//计算矩阵A第x1列y1行开始的列向量模长
		for (i = y1; i < y2; i++) {
			d += A_data[i][x1] * A_data[i][x1];
		}
		d = sqrtf(d);   //现在d是模长。
		//计算向量v。向量v有A.rows-1-y1个元素。
		(beta[0]) = d * (d - A_data[y1][x1]);
		//如果beta等于零，则H为单位矩阵，向量v为零向量。(从beta的计算可知，beta是大于零的)
		if (beta[0] < precision) {
			beta[0] = 1;   	//β不为零是防止下面除法运算出错。 
			for (i = y1; i < y2; i++) {
				v[i - y1] = 0;
			}
		}
		else {				//否则更新Hi。同时更新矩阵A的x1列
			v[0] = A_data[y1][x1] - d;
			A_data[y1][x1] = d;
			for (i = y1 + 1; i < y2; i++) {
				v[i - y1] = A_data[i][x1];
				A_data[i][x1] = 0;
			}
		}

		real_t div_beta = 1 / (beta[0]);
		//计算A左乘Householder矩阵 。A=H(i)*A。从x1+1列开始变换。
		for (j = x1 + 1; j < x2; j++) {	// 
			d = 0;
			for (i = y1; i < y2; i++) {	//计算Householder向量和A列矢的点积
				d = d + v[i - y1] * A_data[i][j];
			}
			d = d * div_beta;
			for (i = y1; i < y2; i++) {
				A_data[i][j] = A_data[i][j] - d * v[i - y1];
			}
		}
	}
	else {
		real_t div_beta = 1 / (beta[0]);
		//计算A左乘Householder矩阵 。A=H(i)*A。从x1列开始变换。
		for (j = x1; j < x2; j++) {
			d = 0;
			for (i = y1; i < y2; i++) {	//计算Householder向量和A列矢的点积
				d = d + v[i - y1] * A_data[i][j];
			}
			d = d * div_beta;
			for (i = y1; i < y2; i++) {
				A_data[i][j] = A_data[i][j] - d * v[i - y1];
			}
		}
	}
}

/// <summary>
/// 豪斯霍尔德变换。
/// 变换后，A = A * H。H是豪斯霍尔德矩阵，以豪斯霍尔德矢量ｖ和参数beta的形式出现。 H=I-(1/beta)v*v^(T)
/// 变换后的A的x1列，y1行以右的元素是0元素，把x1列y1行开始的行向量变成(a,0,0,...)形式的向量。得到的豪斯霍尔德矩阵以向量形式保存在v中。
/// 注意事项：注意数组v的大小，防止越界。
/// 这里是前闭后开。即，前面取得到，后面取不到
/// </summary>
/// <param name="A">矩阵</param>
/// <param name="x1">矩阵豪斯霍尔德变换的起始列</param>
/// <param name="x2">矩阵豪斯霍尔德变换的终止列</param>
/// <param name="y1">矩阵豪斯霍尔德变换的起始行</param>
/// <param name="y2">矩阵豪斯霍尔德变换的终止行</param>
/// <param name="v">豪斯霍尔德矢量，大小是变化列的数量</param>
/// <param name="beta">豪斯霍尔德矢量的参数</param>
/// <param name="isH_ensured">是否已经确认豪斯霍尔德矢量，如果未确认，函数会先计算更新豪斯霍尔德矢量和参数</param>
/// <param name="precision">计算精度</param>
void mat_evd_subf_Householder_right(mat A, int x1, int x2, int y1, int y2, real_t* v, real_t* beta, bool isH_ensured, real_t precision)
{
	int i, j;
	real_t d;
	real_t** A_data = A.data;
	if (x2 > A.cols)x2 = A.cols;
	if (y2 > A.rows)y2 = A.rows;
	if (isH_ensured == false) {
		//先求Householder矩阵的向量
		d = 0;
		//计算矩阵A第x1列y1行开始的行向量模长
		for (i = x1; i < x2; i++) {
			d += A_data[y1][i] * A_data[y1][i];
		}
		d = sqrtf(d);   //现在d是模长。
		//计算向量v。向量v有y2-1-y1个元素。
		beta[0] = d * (d - A_data[y1][x1]);
		//如果beta等于零，则H为单位矩阵，向量v为零向量。(从beta的计算可知，beta是大于零的)
		if (beta[0] < precision) {
			beta[0] = 1;   	//β不为零是防止下面除法运算出错。 
			for (i = x1; i < x2; i++) {
				v[i - x1] = 0;
			}
		}
		else {				//否则更新Hi。同时更新矩阵A的y1行
			v[0] = A_data[y1][x1] - d;
			A_data[y1][x1] = d;
			for (i = x1 + 1; i < x2; i++) {
				v[i - x1] = A_data[y1][i];
				A_data[y1][i] = 0;
			}
		}

		real_t div_beta = 1 / (beta[0]);
		//计算A右乘Householder矩阵 。A=A*H(i)。从y1+1行开始变换。
		for (j = y1 + 1; j < y2; j++) {
			d = 0;
			for (i = x1; i < x2; i++) {			//计算Householder向量和A列矢的点积
				d = d + v[i - x1] * A_data[j][i];
			}
			d = d * div_beta;
			for (i = x1; i < x2; i++) {
				A_data[j][i] = A_data[j][i] - d * v[i - x1];
			}
		}
	}
	else {
		real_t div_beta = 1 / (beta[0]);
		//计算A右乘Householder矩阵 。A=A*H(i)。从y1行开始变换。
		for (j = y1; j < y2; j++) {
			d = 0;
			for (i = x1; i < x2; i++) {			//计算Householder向量和A列矢的点积
				d = d + v[i - x1] * A_data[j][i];
			}
			d = d * div_beta;
			for (i = x1; i < x2; i++) {
				A_data[j][i] = A_data[j][i] - d * v[i - x1];
			}
		}
	}
}

/// <summary>
/// 计算householder矢量
/// </summary>
/// <param name="v">原始矢量</param>
/// <param name="householder">计算得到的householder矢量</param>
/// <param name="beta">householder矢量的beta</param>
/// <param name="vector_length">矢量的长度</param>
/// <param name="precision">计算精度</param>
void mat_get_householder_vector(real_t* v, real_t* householder, real_t* beta, int vector_length, real_t precision)
{
	int i;
	real_t d;
	//先求Householder矩阵的向量
	d = 0;
	//计算矩阵A第0列y1行开始的行向量模长
	for (i = 0; i < vector_length; i++) {
		d += v[i] * v[i];
	}
	d = sqrtf(d);   //现在d是模长。
	//计算向量householder。向量householder有y2-1-y1个元素。
	beta[0] = d * (d - v[0]);
	//如果beta等于零，则H为单位矩阵，向量householder为零向量。(从beta的计算可知，beta是大于零的)
	if (beta[0] < precision) {
		beta[0] = 1;   	//β不为零是防止下面除法运算出错。 
		for (i = 0; i < vector_length; i++) {
			householder[i] = 0;
		}
	}
	else {				//否则更新Hi。同时更新矩阵A的y1行
		householder[0] = v[0] - d;
		for (i = 1; i < vector_length; i++) {
			householder[i] = v[i];
		}
	}
}


/// <summary>
/// 计算非对称矩阵的特征值。这里只计算特征值，不计算特征矢量。
/// </summary>
/// <param name="M">实数方阵</param>
/// <param name="eigen">特征值，特征值的数组要和矩阵的尺寸匹配</param>
/// <param name="precision">精度</param>
/// <returns>迭代次数</returns>
int mat_EVD_only_Eigenvalue(mat M, complex* eigen, real_t precision)
{
	mat D; D.rows = M.rows; D.cols = M.rows;

	mat_allocate(&D, M.rows, M.rows);
	vec sub_diag;
	vec_allocate(&sub_diag, M.rows);

	real_t v[3];
	real_t householder[3];
	real_t beta[1];

	int times = -1;
	if (EVD_Hessenberg_Householder_without_H(M, D, precision)) {

		//现在X是单位正交矩阵，D是上hessenberg矩阵。
		//使用双重步位移
	//	if (precision < 0.0001f) precision = 0.0001f;  //精度不能太小。

		real_t** D_data = D.data;
		int i = 0;
//		int j = 0;
		int sm = 0, sn = 0; //记录qr分解的子矩阵的位置。从sm行列到sn行列。
		bool is_diag = false;
		bool flag = false;
		bool isget_mn = false;
		int sub_times;  //记录子循环的次数。

		// 通过计算次对角元的变化来判断是否收敛。如果连续两次迭代，次对角元素最大的变化小于0.001，则认为已经收敛了。

		// 记录次对角元素。
		for (i = 0; i < M.cols - 1; i++) {
			sub_diag.data[i] = D.data[i + 1][i];
		}
		real_t max_err0 = 1000;
//		real_t max_err1 = 1000;
//		real_t max_err2 = 1000;

		int Th = powf(M.rows, 1.5f);
		while (is_diag == false) //下面分块进行QR分解迭代。
		{
			//设计思路:设初始flag为0，依次检测斜对角元素，为零flag置1，否则置0。检测上升沿和下降沿。
			//上升沿设置sm，下降沿设置sn。[sm，sn]是QR迭代的子三对角矩阵区域。设置了mn即可开始迭代。sm，sn两端都是可以取的。
			//如果无下降沿，则sn=m.cols-1。如无上升沿，则斜对角线都是零元，可结束迭代。

			// 如果迭代次数太大，需要放宽精度。
			if (Th > M.rows) {
				if (times % (Th + 10) == Th) {
					precision *= 2;
				}
			}
			//初始化标志位
			flag = false;
			isget_mn = false;
			sub_times = 0;
			for (i = 1; i < M.cols; i++) {
				real_t ep = Max(precision, precision * (ABS(D_data[i][i]) + ABS(D_data[i - 1][i - 1])));
				if (ABS(D_data[i][i - 1]) > ep) {  //非零
					if (flag == false) {	//检测上升沿。
						sm = i - 1;
						flag = true;
					}
					if (flag == true && i == M.cols - 1) {
						sn = i;
						isget_mn = true;
						flag = false;
					}
				}
				else {
					if (flag == true) {		//检测下降沿
						sn = i - 1;
						isget_mn = true;
					}
					flag = false;
				}
				if (isget_mn && sn - sm > 1) {  // 出现不是2*2的 对角矩阵元

					//QR分解迭代，双重步位移法
					int m, n = sn;
					{
						m = n - 1;
						real_t s, t;
						s = D.data[m][m] + D.data[n][n];
						t = D.data[m][m] * D.data[n][n] - D.data[m][n] * D.data[n][m];
						v[0] = D.data[sm][sm] * D.data[sm][sm] + D.data[sm][sm + 1] * D.data[sm + 1][sm] - s * D.data[sm][sm] + t;
						v[1] = D.data[sm + 1][sm] * (D.data[sm][sm] + D.data[sm + 1][sm + 1] - s);
						v[2] = D.data[sm + 1][sm] * D.data[sm + 2][sm + 1];

						int k;
						for (k = sm; k <= sn - 2; k++)
						{
							mat_get_householder_vector(v, householder, beta, 3, precision);
							int q = Max(sm, k - 1);  //Householder_left 起始列数  ，k是起始行数。
							int r = Min(k + 3, sn);  // 终止行数 (取)
							mat_evd_subf_Householder_left(D, q, D.cols, k, k + 3, householder, beta, true, precision);
							mat_evd_subf_Householder_right(D, k, k + 3, 0, r + 1, householder, beta, true, precision);
							v[0] = D.data[k + 1][k];
							v[1] = D.data[k + 2][k];
							if (k <= sn - 3)
								v[2] = D.data[k + 3][k];
						}
						mat_get_householder_vector(v, householder, beta, 2, precision);
						mat_evd_subf_Householder_left(D, sn - 2, sn + 1, sn - 1, sn + 1, householder, beta, true, precision);
						mat_evd_subf_Householder_right(D, sn - 1, sn + 1, 0, sn + 1, householder, beta, true, precision);

					}
					//迭代一次后就重置标记
					isget_mn = false;
					sub_times++;
					//	printfs(" sub_times =  " << sub_times << "  。 \r\n";
				}
			}
			//如果扫描了所有斜对角线元素，都没有确定一个mn，则说明全是零元，无上升沿。分解成对角矩阵了。
			if (sub_times == 0) {
				is_diag = true;
			}
			times++;

			//计算最大的误差值。
//			max_err2 = max_err1;
//			max_err1 = max_err0;
			max_err0 = 0;
			for (i = 0; i < M.cols - 1; i++) {
				max_err0 = Max(max_err0, ABS(ABS(sub_diag.data[i]) - ABS(D.data[i + 1][i])));
			}
			// 调试发现，max_err012 都在0.5左右徘徊，有一次会降低到0.0004，但随后又升到0.5…… 后面应该是舍入误差了…… 还是防止精度过小，比较有效

			// 记录次对角元素。
			for (i = 0; i < M.cols - 1; i++) {
				sub_diag.data[i] = D.data[i + 1][i];
			}
				//printfs("    this is    %d  times   \r\n",times);
				//mat_print(D);
		}

		// 下面计算特征值。
		for (i = 0; i < M.cols; i++) {
			if (i < M.cols - 1) {  //检测是否 2*2 矩阵。
				real_t ep = Max(precision, precision * (ABS(D_data[i + 1][i + 1]) + ABS(D_data[i][i])));
				//	if (times > D.cols * D.cols * 2) ep = (real_t)(times) * (D.cols * D.cols * 2); // 如果每行每列至少迭代了4次，还在迭代，得放宽收敛条件了，不然无法退出循环。

			//	bool is_larger_then_0 = false;
				real_t a, b, c;

				a = (i > 0) ? ABS(D_data[i][i - 1]) : 0;
				b = ABS(D_data[i + 1][i]);
				c = (i < M.cols - 2) ? ABS(D_data[i + 2][i + 1]) : 0;

				if (b > a && b > c) { //非零
					real_t r1 = D_data[i][i] + D_data[i + 1][i + 1];
					real_t r2 = D_data[i][i] * D_data[i + 1][i + 1] - D_data[i][i + 1] * D_data[i + 1][i];

					real_t delta = r1 * r1 - 4 * r2;
					if (delta < 0) {
						delta = sqrtf(-delta);
						eigen[i].imag = 0.5 * delta;
						eigen[i].real = 0.5 * (r1);
						i++;
						eigen[i].imag = -0.5 * delta;
						eigen[i].real = 0.5 * (r1);
					}
					else {
						delta = sqrtf(delta);
						eigen[i].imag = 0;
						eigen[i].real = 0.5 * (r1 + delta);
						i++;
						eigen[i].imag = 0;
						eigen[i].real = 0.5 * (r1 - delta);
					}
				}
				else {
					eigen[i].imag = 0;
					eigen[i].real = D_data[i][i];
				}
			}
			else {
				eigen[i].imag = 0;
				eigen[i].real = D_data[i][i];
			}
		}

	}

	// mat_print(D);
	// printfs(" evd finish  \r\n \r\n ");
	mat_deallocate(&D);
	vec_deallocate(&sub_diag);
	return times;
}
//

/// <summary>
/// 已知 a_0 + a_1*x + a_2*x^2 + ... + a_(n-1)*x^(n-1) + x^n = 0
/// 测试结果：和matlab的roots函数相比，大概有0.01的误差。勉强能用。  matlab的矩阵计算除以零元，及对精度的处理有优化，可以减小舍入误差。
/// 而我这里水平和条件有限，先做到这样吧。
/// real_t是双精度浮点型的时候，计算的精度比较高，用解出来的根代进原多项式，也比较接近0。
/// 如果real_t是单精度类型，把根代进多项式，当根的模长比较大的时候，值会比较大（精度不够）。高阶多项式偏离零点比较远的时候，会很不稳定。
/// </summary>
/// <param name="coefficient"> 系数 a0,a1,a2……的数组，数组长度是rank</param>
/// <param name="root">根,即函数x的解，数组长度是rank</param>
/// <param name="rank">方程的阶数。 最高阶的系数默认是1</param>
/// <returns>迭代次数</returns>
int roots(real_t* coefficient, complex* root, int rank, real_t precision)
{
	mat A;
	mat_allocate(&A, rank, rank);
	mat_init(A, 0);
	int i;
	for (i = 0; i < rank; i++) {
		A.data[0][i] = -coefficient[rank - i - 1];
	}
	for (i = 1; i < rank; i++) {
		A.data[i][i - 1] = 1;
	}
	int a = mat_EVD_only_Eigenvalue(A, root, precision);
	mat_deallocate(&A);
	return a;
}


///////////////   测试         ////////////////////////////////////


/// 测试向量相关的运算
void testVector(void)
{
	int Length = 8;

	vec A;
	vec B;
	vec C;
	vec D;
	vec E;
	vec_allocate(&A, Length);
	vec_allocate(&B, Length);
	vec_allocate(&C, Length);
	vec_allocate(&D, Length);
	vec_allocate(&E, Length);

	vec_randomize(A, 314159,0,1);
	vec_randomize(B, 31415,0,1);

	vec_add(A, B, C);
	vec_sub(C, B, D);
	vec_sub(D, A, E);


	if (vec_isZero(E, Epsilon) && vec_isEqual(A, D, Epsilon)) {
		printfs("vec_add success \n   \r\n");
	}
	else {
		printfs("vec_add  failed \n   \r\n");
	}

	printfs("  A mod square %f \r\n ", vec_length(A) * vec_length(A));
	//验证乘法。A的模长*2 = (2*A)点乘A
	vec_mul(A, 2, B);
	if (ABS(vec_dotProduct(A, B) - 2 * vec_length(A) * vec_length(A)) < Epsilon) {
		printfs("vec_dotProduct success \n   \r\n");
	}
	else {
		printfs("vec_dotProduct failed \n   \r\n");
	}

	vec_deallocate(&A);
	vec_deallocate(&B);
	vec_deallocate(&C);
	vec_deallocate(&D);
	vec_deallocate(&E);

}


//用于验证矩阵类加减乘，求逆，QR分解，奇异值，特征值分解的可靠性。利用逆运算来进行验证。
void test_Matrix_add_sub(void)
{
	int row = 6, col = 6;

	mat m1;
	mat m2;
	mat m3;
	mat m4;
	mat m5;
	mat_allocate(&m1, row, col);
	mat_allocate(&m2, row, col);
	mat_allocate(&m3, row, col);
	mat_allocate(&m4, row, col);
	mat_allocate(&m5, row, col);

	//验证加减法
	mat_randomize(m1, 314159,0,1);
	mat_randomize(m2, 31415,0,1);
	mat_add(m1, m2, m3);
	mat_sub(m3, m2, m4);
	if (mat_isEqual(m1, m4, Epsilon)) {
		printfs( "mat_add success \r\n");
	}
	else {
		printfs( "mat_add failed  \r\n");
	}

	//验证乘除法
	mat_randomize(m1, 314159, 0, 1);
	mat_randomize(m2, 31415, 0, 1);
	mat_mul(m1, m2, m3);
	mat_inversion(m2, m4, Epsilon);
	mat_mul(m3, m4, m5);
	if (mat_isEqual(m1, m5, Epsilon * 100)) {
		printfs("mat_inversion success \r\n");
	}
	else {
		printfs("mat_inversion failed  \r\n");
	}

	mat_deallocate(&m1);
	mat_deallocate(&m2);
	mat_deallocate(&m3);
	mat_deallocate(&m4);
	mat_deallocate(&m5);

}

bool test_Matrix_mul_div(void)
{
	return 1;
}

bool test_Matrix_det(void)
{
	return 1;
}

void test_Matrix_QR_Givens(void)
{
	int r[3] = { 11,13,17 };
	int c[3] = { 17,13,11 };

	printfs(" \r\n  Givens QR test \r\n");
	for (int i = 0; i < 3; i++) {
		int row = r[i];
		int col = c[i];
		mat m;
		mat Q;
		mat R;
		mat QR;
		mat QTQ;
		mat QT;

		mat_allocate(&m, row, col);
		mat_allocate(&Q, row, row);
		mat_allocate(&R, row, col);
		mat_allocate(&QR, row, col);
		mat_allocate(&QTQ, row, row);
		mat_allocate(&QT, row, row);

		mat_randomize(m, 314159, -2, 2);
		{
			//首尾行置零+中间行
			for (int j = 0; j < col; j++) {
				m.data[0][j] = 0;
				m.data[row / 2][j] = 0;
				m.data[row / 2 + 1][j] = 0;
				m.data[row - 1][j] = 0;
			}
			//首尾列置零+中间列
			for (int i = 0; i < row; i++) {
				m.data[i][0] = 0;
				m.data[i][col / 2] = 0;
				m.data[i][col / 2 + 1] = 0;
				m.data[i][col - 1] = 0;
			}
		}
		mat_QR_Givens(m, Q, R, Epsilon);

		mat_mul(Q, R, QR);
		if (mat_isEqual(m, QR,Epsilon)) {
			printfs( "QR success \r\n");
		}
		else {
			printfs( "QR failed  \r\n");
			mat_print(m);
			mat_print(Q);
			mat_print(R);
			mat_print(QR);
		}
		mat_transpose_to(Q, QT);
		mat_mul(QT, Q, QTQ);
		if (mat_isEye(QTQ,Epsilon)) {
			printfs( "Q mul QT is eye,success \r\n");
		}
		else {
			printfs( "Q mul QT is not eye,failed  \r\n");
			mat_print(QTQ);
		}
		if (mat_isUpperTriangle(R,Epsilon)) {
			printfs( "R isUpperTriangle , success  \r\n");
		}
		else {
			printfs( "R is not UpperTriangle , failed     \r\n");
			mat_print(R);
		}

		mat_deallocate(&m);
		mat_deallocate(&Q);
		mat_deallocate(&R);
		mat_deallocate(&QR);
		mat_deallocate(&QTQ);
		mat_deallocate(&QT);
	}
}
void test_Matrix_QR_Schmidt(void)
{
	int r[3] = { 11,13,17 };
	int c[3] = { 17,13,11 };
	//int r[3] = { 11,12,13 };
	//int c[3] = { 13,12,11 };

	printfs(" \r\n  Schmidt QR test \r\n");
	for (int i = 0; i < 3; i++) {
		int row = r[i];
		int col = c[i];
		mat m;
		mat Q;
		mat R;
		mat QR;
		mat QTQ;
		mat QT;

		mat_allocate(&m, row, col);// 格拉姆施密特过程不要求原矩阵是方阵，只要Q的行列数和原矩阵m一致，R是与列数相同的方阵。
		mat_allocate(&Q, row, col);
		mat_allocate(&R, col, col);
		mat_allocate(&QR, row, col);
		mat_allocate(&QTQ, col, col);
		mat_allocate(&QT, col, row);

		mat_randomize(m, 34159, -2, 2);
		if(1)
		{
			//首尾行置零+中间行
			for (int j = 0; j < col; j++) {
				m.data[0][j] = 0;
				m.data[row / 2][j] = 0;
				m.data[row / 2 + 1][j] = 0;
				m.data[row - 1][j] = 0;
			}
			//首尾列置零+中间列
			for (int i = 0; i < row; i++) {
				m.data[i][0] = 0;
				m.data[i][col / 2] = 0;
				m.data[i][col / 2 + 1] = 0;
				m.data[i][col - 1] = 0;
			}
		}
		mat_QR_Schmidt(m, Q, R, Epsilon);  

		mat_mul(Q, R, QR);
		if (mat_isEqual(m, QR, Epsilon)) {
			printfs("QR success \r\n");
		}
		else {
			printfs("QR failed  \r\n");
			mat_print(m);
			mat_print(Q);
			mat_print(R);
			mat_print(QR);
		}
		mat_transpose_to(Q, QT);
		mat_mul(QT, Q, QTQ);
		if (mat_isDiagonal_01(QTQ, Epsilon)) {
			printfs("Q mul QT is eye,success \r\n");
		}
		else {
			printfs("Q mul QT is not eye,failed  \r\n");

			printfs("Q   \r\n");
			mat_print(Q);
			printfs("R   \r\n");
			mat_print(R);
			printfs(" QT  \r\n");
			mat_print(QT);
			printfs("Q mul QT  \r\n");
			mat_print(QTQ);

		}
		if (mat_isUpperTriangle(R, Epsilon)) {
			printfs("R isUpperTriangle , success  \r\n");
		}
		else {
			printfs("R is not UpperTriangle , failed     \r\n");
			mat_print(R);
		}

		mat_deallocate(&m);
		mat_deallocate(&Q);
		mat_deallocate(&R);
		mat_deallocate(&QR);
		mat_deallocate(&QTQ);
		mat_deallocate(&QT);
	}
}
void test_Matrix_QR_Householder(void)
{
	int r[3] = { 11,13,17 };
	int c[3] = { 17,13,11 };

	printfs(" \r\n  Householder QR test \r\n");
	for (int i = 0; i < 3; i++) {
		int row = r[i];
		int col = c[i];
		mat m;
		mat Q;
		mat R;
		mat QR;
		mat QTQ;
		mat QT;

		mat_allocate(&m, row, col);  //不要求m矩阵是方阵。A尺寸与m矩阵一样。H是和行数一样的方阵。
		mat_allocate(&Q, row, row);
		mat_allocate(&R, row, col);
		mat_allocate(&QR, row, col);
		mat_allocate(&QTQ, row, row);
		mat_allocate(&QT, row, row);


		mat_randomize(m, 314159,-2,2);
		{
			//首尾行置零+中间行
			int setvalue = 1;

			for (int j = 0; j < col; j++) {
				m.data[0][j] = setvalue;
				m.data[row / 2][j] = setvalue;
				m.data[row / 2 + 1][j] = setvalue;
				m.data[row - 1][j] = setvalue;
			}
			//首尾列置零+中间列
			for (int i = 0; i < row; i++) {
				m.data[i][0] = setvalue;
				m.data[i][col / 2] = setvalue;
				m.data[i][col / 2 + 1] = setvalue;
				m.data[i][col - 1] = setvalue;
			}
		}


		mat_QR_Householder(m, Q, R, Epsilon);

		mat_mul(Q, R, QR);
		if (mat_isEqual(m, QR, Epsilon)) {
			printfs("QR success \r\n");

			//printfs("m  \r\n");
			//mat_print(m);
			//printfs("Q  \r\n");
			//mat_print(Q);
			//printfs("R  \r\n");
			//mat_print(R);
			//printfs("QR  \r\n");
			//mat_print(QR);

		}
		else {
			printfs("QR failed  \r\n");

			printfs("m  \r\n");
			mat_print(m);
			printfs("Q  \r\n");
			mat_print(Q);
			printfs("R  \r\n");
			mat_print(R);
			printfs("QR  \r\n");
			mat_print(QR);
		}
		mat_transpose_to(Q, QT);
		mat_mul(QT, Q, QTQ);
		if (mat_isEye(QTQ, Epsilon)) {
			printfs("Q mul QT is eye,success \r\n");
		}
		else {
			printfs("Q mul QT is not eye,failed  \r\n");
			mat_print(QTQ);
		}
		if (mat_isUpperTriangle(R, Epsilon)) {
			printfs("R isUpperTriangle , success  \r\n");
		}
		else {
			printfs("R is not UpperTriangle , failed     \r\n");
			mat_print(R);
		}

		mat_deallocate(&m);
		mat_deallocate(&Q);
		mat_deallocate(&R);
		mat_deallocate(&QR);
		mat_deallocate(&QTQ);
		mat_deallocate(&QT);
	}
}
//
void test_Matrix_LinearLeastSquare(void)
{
	//m*x=v;已知m与v，求x。
	//验证最小二乘的方法:m*x是最接近v的。所以v-m*x是垂直于m*x的。  v-m*x 也是垂直于所有 m 的列向量的。
	int r[3] = { 11,13,17 };
	int c[3] = { 17,13,11 };
	//int r[3] = { 11,12,13 };
	//int c[3] = { 13,12,11 };

	printfs(" \r\n  Schmidt QR test \r\n");
	for (int i = 0; i < 3; i++) {
		int row = r[i];
		int col = c[i];

		mat m1;
		vec v; 
		vec x; 
		vec mx;
		vec v_mx; 
		vec m1T_mul_v_mx;

		//m*x=v;已知m与v，求x。
		mat_allocate(&m1, row, col);
		vec_allocate(&x, col);
		vec_allocate(&v, row);
		vec_allocate(&mx, row);
		vec_allocate(&v_mx, row);
		vec_allocate(&m1T_mul_v_mx, col);

		mat_randomize(m1, 314159,-2,2);
		vec_randomize(v, 27182,-2,2);
		{
			//首尾行置零+中间行
			for (int j = 0; j < col; j++) {
				m1.data[0][j] = 0;
				m1.data[row / 2][j] = 0;
				m1.data[row / 2 + 1][j] = 0;
				m1.data[row - 1][j] = 0;
			}
			//首尾列置零+中间列
			for (int i = 0; i < row; i++) {
				m1.data[i][0] = 0;
				m1.data[i][col / 2] = 0;
				m1.data[i][col / 2 + 1] = 0;
				m1.data[i][col - 1] = 0;
			}
		}

		mat_LinearEquation_LinearLeastSquare(m1, x, v, true, Epsilon);
		mat_mul_vec(m1, x, mx);
		vec_sub(v, mx, v_mx);
		if (ABS(vec_dotProduct(v_mx, mx)) < Epsilon * mx.size ) {
			printfs("\r\n  mat_LinearEquation_LinearLeastSquare success ");
		}
		else {
			matT_mul_vec(m1, v_mx, m1T_mul_v_mx);

			printfs("\r\n  mat_LinearEquation_LinearLeastSquare failed ");
			vec_print(m1T_mul_v_mx);
			printfs("\r\n %f   ====  ", vec_dotProduct(v_mx, mx));
			vec_print(x);
			vec_print(mx);
			vec_print(v_mx);

		}

		mat_deallocate(&m1);
		vec_deallocate(&v);
		vec_deallocate(&x);
		vec_deallocate(&mx);
		vec_deallocate(&v_mx);
		vec_deallocate(&m1T_mul_v_mx);


	}
}

//
void test_Matrix_EVD(void) 
{
	int row = 114;
	int col = 114;
	mat m;
	mat Q;
	mat R;
	mat QR;
	mat QTQ;
	mat QT;

	mat_allocate(&m, row, col);  
	mat_allocate(&Q, row, row);
	mat_allocate(&R, row, col);
	mat_allocate(&QR, row, col);
	mat_allocate(&QTQ, row, row);
	mat_allocate(&QT, row, row);

	mat_randomize(m, 314159, -2, 2);
	for (int i = 0; i < row; i += 1) {
		for (int j = 0; j < col; j += 1) {
			m.data[j][i] = m.data[i][j];
		}
	}
	//中间行列置零
	if (1) {
		for (int i = 0; i < row; i += 1) {
			m.data[i][col / 2] = 0;
			m.data[col / 2][i] = m.data[i][col / 2];

			m.data[i][col / 2 + 1] = 0;
			m.data[col / 2 + 1][i] = m.data[i][col / 2 + 1];
		}
	}

	mat_EVD_Symmetric(m, Q, R, Epsilon);

	mat_mul(Q, R, QR);
	mat_transpose_to(Q, QT);
	mat_mul(Q, R, QR);
	mat_mul(QR, QT, QTQ);

	if (mat_isEqual(m, QTQ, Epsilon * 10)) {
		printfs( "EVD success \r\n");
	}
	else { 
		printfs( "EVD faled \r\n");
		mat_print(m);
		mat_print(Q);
		mat_print(R);
		mat_print(QTQ);

	}

	mat_transpose_to(Q, QT);
	mat_mul(QT, Q, QTQ);
	if (mat_isEye(QTQ, Epsilon * 10)) {
		printfs( "QTQ isEye ,success   \r\n");
	}
	else {
		printfs( "QTQ is not Eye ,failed \r\n");
		mat_print(QTQ);
	}
	if (mat_isLowerTriangle(R, Epsilon * 10) && mat_isUpperTriangle(R, Epsilon * 10)) {
		printfs( "R is diagonal ,success \r\n");
	}
	else {
		printfs( "R is not diagonal ,failed \r\n");
		mat_print(R);
	}


	mat_deallocate(&m);  //不要求m矩阵是方阵。A尺寸与m矩阵一样。H是和行数一样的方阵。
	mat_deallocate(&Q);
	mat_deallocate(&R);
	mat_deallocate(&QR);
	mat_deallocate(&QTQ);
	mat_deallocate(&QT);

}

//
void test_Matrix_SVD(void)
{
	int r[3] = { 111,113,117 };
	int c[3] = { 117,113,111 };

	//int r[3] = { 13,6,7 };
	//int c[3] = { 7,6,13 };

	printfs(" \r\n  test_Matrix_SVD \r\n");
	for (int i = 0; i < 3; i++) {


		int row = r[i], col = c[i];
		int times = 0;
		/*
			real_t data__m [row][col];
			real_t data__cm[row][col];
			real_t data__s[row][col];

			real_t data__u[row][row];
			real_t data__vt[col][col];
			real_t data__ut[row][row];
			real_t data__v[col][col];
			real_t data__utu[row][row];
			real_t data__vtv[col][col];
		*/
		mat m, cm,cm2, s, u, ut, utu, v, vt, vtv;

		mat_allocate(&m, row, col);
		mat_allocate(&cm, row, col);
		mat_allocate(&cm2, row, col);
		mat_allocate(&s, row, col);
		mat_allocate(&u, row, row);
		mat_allocate(&ut, row, row);
		mat_allocate(&utu, row, row);
		mat_allocate(&v, col, col);
		mat_allocate(&vt, col, col);
		mat_allocate(&vtv, col, col);

		mat_randomize(m, 314159, 0, 1);
		if (1)
		{
			//首尾行置零+中间行
			for (int j = 0; j < col; j++) {
				m.data[0][j] = 0;
				m.data[row / 2][j] = 0;
				m.data[row / 2 + 1][j] = 0;
				m.data[row - 1][j] = 0;
			}
			//首尾列置零+中间列
			for (int i = 0; i < row; i++) {
				m.data[i][0] = 0;
				m.data[i][col / 2] = 0;
				m.data[i][col / 2 + 1] = 0;
				m.data[i][col - 1] = 0;
			}
		}


		//mat_copyto(m, s);
		//times = SVD_subfunction1(m, u, s, vt, Epsilon);
		//SVD_subfunction3(u, s, Epsilon);
		//SVD_subfunction2(u, s, vt, 0, 5, Epsilon);

		 times = mat_SVD(m, u, s, vt, Epsilon);


		printfs(" \r\n %d  \r\n", times);
	//	mat_print(s);

		mat_mul(u, s, cm);
		mat_mul(cm, vt, cm2);
		if (mat_isEqual(m, cm2, Epsilon*10000)) {
			printfs(" \r\n ============== SVD                     success \r\n");

			//printfs(" m \r\n");
			//mat_print(m);
			//printfs(" cm2 \r\n");
			//mat_print(cm2);
			//printfs(" u \r\n");
			//mat_print(u);
			//printfs(" s \r\n");
			//mat_print(s);
			//printfs(" vt \r\n");
			//mat_print(vt);

		}
		else {
			printfs(" \r\n ========== SVD                      failed \r\n");

			printfs(" m \r\n");
			mat_print(m);
			printfs(" cm2 \r\n");
			mat_print(cm2);
			printfs(" u \r\n");
			mat_print(u);
			printfs(" s \r\n");
			mat_print(s);
			printfs(" vt \r\n");
			mat_print(vt);
		}

		mat_transpose_to(u, ut);
		mat_mul(ut, u, utu);
		if (mat_isEye(utu, Epsilon)) {
			printfs("u is eye ,success \r\n");
		}
		else {
			printfs("u is not eye ,failed \r\n");
			mat_print(utu);
		}

		mat_transpose_to(vt, v);
		mat_mul(vt, v, vtv);
		if (mat_isEye(vtv, Epsilon)) {
			printfs("v is eye,success \r\n");
		}
		else {
			printfs("v is not eye,failed \r\n");
			mat_print(vtv);
		}


		mat_deallocate(&m);
		mat_deallocate(&cm);
		mat_deallocate(&cm2);
		mat_deallocate(&s);
		mat_deallocate(&u);
		mat_deallocate(&ut);
		mat_deallocate(&utu);
		mat_deallocate(&v);
		mat_deallocate(&vt);
		mat_deallocate(&vtv);

	}
}

void testMatrix(void)
{
	test_Matrix_add_sub();
	test_Matrix_mul_div();
	test_Matrix_det();

	test_Matrix_QR_Givens();
	test_Matrix_QR_Schmidt();
	test_Matrix_QR_Householder();
	test_Matrix_LinearLeastSquare();
	test_Matrix_EVD();

	test_Matrix_SVD();

}

//


#define LenR 13
void test_roots()  
{
	printfs(" \r\n test roots \r\n ");
//	real_t a[LenR] = { 1,2,3,4,5,6,7,8,9 };

//	real_t a[LenR] = { 9,8,7,6,5,4,3,2,1 };

//	real_t a[LenR] = { 5,3,4,5,6,7,8,2,16 };

	real_t a[LenR] = { 2.5,2.9,6.8,40.2,6 };

	complex c[LenR];
	int times = roots(a, c, LenR, Epsilon);

	printfs("\r\n times : %d", times);

	for (int i = 0; i < LenR; i++) {
		printfs(" \r\n root :");
		complex_print(c[i]);
	}

	// 计算函数值，反向验证

	for (int k = 0; k < LenR; k++) {
		complex x = c[k];

		complex y = { a[0],0 };
		complex temp;
		complex xn = x;
		for (int i = 1; i < LenR; i++) {
			temp = complex_mul_real(xn, a[i]);
			y = complex_add(y, temp);
			xn = complex_mul(x, xn);
		}
		y = complex_add(y, xn);

		printfs(" \r\n root result :");
		complex_print(y);


	}

}



void test_mat_house()
{
//#define row 5
//#define col 5

	const int row = 5, col = 5;
	mat A;
	mat_allocate(&A, row, col);
	mat_randomize(A, 314159, -2, 2);
	for (int i = 0; i < row; i += 1) {
		for (int j = 0; j < col; j += 1) {
			A.data[j][i] = A.data[i][j];
		}
	}

	real_t v[5];
	real_t householder[5];
	real_t beta[1];
	real_t precision = Epsilon;
	
	for (int i = 1; i < row; i++) {

		for (int j = i; j < row; j++) {
			v[j - i] = A.data[j][i - 1];
		}

		mat_get_householder_vector(v, householder, beta, row - i, precision);
		mat_evd_subf_Householder_left(A, i - 1, col, i, row, householder, beta, true, precision);
		mat_evd_subf_Householder_right(A, i, row, i - 1, col, householder, beta, true, precision);
		mat_print(A);

	}

	printfs("   \r\n  test_mat_house finish  ");

	mat_deallocate(&A);


}


void test_mat_EVD_only_Eigenvalue(void)
{
	test_roots();
	//  return;

	int row = 23;
	int col = 23;
	complex eigen[100];
	mat m;
	mat m1;
	mat Q;
	mat R;
	mat QR;
	mat QTQ;
	mat QT;

	mat_allocate(&m, row, col);
	mat_allocate(&m1, row, col);
	mat_allocate(&Q, row, row);
	mat_allocate(&R, row, col);
	mat_allocate(&QR, row, col);
	mat_allocate(&QTQ, row, row);
	mat_allocate(&QT, row, row);

	mat_randomize(m, 314159, -2, 2);

	for (int i = 0; i < row; i += 1) {
		for (int j = 0; j < col; j += 1) {
			m.data[j][i] = m.data[i][j];
		}
	}
	
	// 中间行列置零
	if (1) {
		for (int i = 0; i < row; i += 1) {
			m.data[i][col / 2] = 0;
			m.data[col / 2][i] = m.data[i][col / 2];
			m.data[i][col / 2 + 1] = 0;
			m.data[col / 2 + 1][i] = m.data[i][col / 2 + 1];
		}
	}
	int times = 0;

	times = mat_EVD_Symmetric(m, Q, R, Epsilon);
	// mat_print(R);
	mat_copyfrom(m1, m);

	times = mat_EVD_only_Eigenvalue(m1, eigen, Epsilon);   //    这个很灵活多变……  
//	times = mat_EVD_only_Eigenvalue(m1, eigen, 0.001);   //    这个很灵活多变……  
	printfs("\r\n times : %d", times);
	for (int i = 0; i < col; i++) {
		printfs(" \r\n egin :");
		complex_print(eigen[i]);
	}
//	mat_EVD_only_Eigenvalue(m1, Epsilon*100);
//	mat_print(m1);

	mat_mul(Q, R, QR);
	mat_transpose_to(Q, QT);
	mat_mul(Q, R, QR);
	mat_mul(QR, QT, QTQ);


	if (mat_isEqual(m, QTQ, Epsilon * 10)) {
		printfs("EVD success \r\n");
		mat_print(m);
		mat_print(Q);
		mat_print(R);
		mat_print(QTQ);
	}
	else {
		printfs("EVD faled \r\n");
		//mat_print(m);
		//mat_print(Q);
		//mat_print(R);
		//mat_print(QTQ);

	}

	mat_transpose_to(Q, QT);
	mat_mul(QT, Q, QTQ);
	if (mat_isEye(QTQ, Epsilon * 10)) {
		printfs("QTQ isEye ,success   \r\n");
	}
	else {
		printfs("QTQ is not Eye ,failed \r\n");
	//	mat_print(QTQ);
	}
	if (mat_isLowerTriangle(R, Epsilon * 10) && mat_isUpperTriangle(R, Epsilon * 10)) {
		printfs("R is diagonal ,success \r\n");
	}
	else {
		printfs("R is not diagonal ,failed \r\n");
	//	mat_print(R);
	}


	mat_deallocate(&m);  //
	mat_deallocate(&m1);  
	mat_deallocate(&Q);
	mat_deallocate(&R);
	mat_deallocate(&QR);
	mat_deallocate(&QTQ);
	mat_deallocate(&QT);

}




